# -*- coding: utf-8 -*-
"""Full-local backend: GRD from CDSE / ASF / Planetary Computer; openEO-aligned σ⁰ PWTT."""

import math
import os
import threading
from datetime import datetime
import numpy as np
from typing import Optional, Tuple

from qgis.core import QgsSettings

from .base_backend import PWTTBackend
from .downloader import get_token, search_s1_grd, download_product, find_vv_vh_in_safe
from .local_numpy_ops import (
    convolve2d_edge,
    openeo_style_p_value_bound,
    welford_init,
    welford_sample_variance,
    welford_update,
)
from .utils import wkt_to_bbox

LOCAL_SOURCE_CDSE = "cdse"
LOCAL_SOURCE_ASF = "asf"
LOCAL_SOURCE_PC = "pc"

# Default cap on scenes per period (pre/post). Higher values track openEO (all acquisitions
# in the window) better; can override via QgsSettings PWTT/local_max_scenes_per_period.
_DEFAULT_MAX_SCENES_PER_PERIOD = 24


def _max_scenes_per_period() -> int:
    s = QgsSettings()
    s.beginGroup("PWTT")
    v = s.value("local_max_scenes_per_period", _DEFAULT_MAX_SCENES_PER_PERIOD)
    s.endGroup()
    try:
        n = int(v)
    except (TypeError, ValueError):
        n = _DEFAULT_MAX_SCENES_PER_PERIOD
    return max(1, min(n, 80))


def _parse_relative_orbit_from_name(name: str) -> Optional[int]:
    """Derive Sentinel-1 relative orbit from a SAFE product name.

    SAFE filename form: S1A_IW_GRDH_1SDV_20220301T041201_..._045512_057123_....
    Field 7 (0-indexed 6) is the absolute orbit. S1A: rel = ((abs-73)%175)+1;
    S1B: rel = ((abs-27)%175)+1.
    """
    if not name:
        return None
    parts = name.split("_")
    if len(parts) < 7:
        return None
    mission = parts[0]
    try:
        abs_orbit = int(parts[6])
    except (TypeError, ValueError):
        return None
    if mission == "S1A":
        return ((abs_orbit - 73) % 175) + 1
    if mission == "S1B":
        return ((abs_orbit - 27) % 175) + 1
    return None


# Default ground-range pixel spacing for AOI warp (~Sentinel-1 GRD IW).
_AOI_RESOLUTION_M = 10.0
# Expand AOI in projected metres so edge pixels are not clipped by reprojection.
_AOI_MARGIN_M = 40.0

# Keep progress lines readable in the Jobs dock (full text still in exceptions / job log).
_PROGRESS_ERR_MAX_LEN = 200


def _is_identity_pixel_transform(t) -> bool:
    """True when *t* is the default 1×1 pixel grid (no real geotransform)."""
    return (
        abs(t.a - 1.0) < 1e-9
        and abs(t.e - 1.0) < 1e-9
        and abs(t.b) < 1e-12
        and abs(t.d) < 1e-12
        and abs(t.c) < 1e-12
        and abs(t.f) < 1e-12
    )


def _aoi_utm_grid(
    west: float,
    south: float,
    east: float,
    north: float,
    resolution_m: float = _AOI_RESOLUTION_M,
    margin_m: float = _AOI_MARGIN_M,
) -> Tuple[object, object, int, int]:
    """Build a UTM raster grid covering the WGS84 bounding box.

    Returns ``(dst_crs, dst_transform, height, width)`` suitable for warping all scenes
    to a common pixel grid (avoids full-frame loads and fixes GCP-only products).
    """
    from rasterio.crs import CRS
    from rasterio.transform import from_bounds
    from rasterio.warp import transform_bounds

    lon_c = (west + east) / 2.0
    lat_c = (south + north) / 2.0
    zone = int((lon_c + 180.0) // 6) + 1
    epsg = (32600 if lat_c >= 0 else 32700) + zone
    dst_crs = CRS.from_epsg(epsg)
    left, bottom, right, top = transform_bounds(
        "EPSG:4326", dst_crs, west, south, east, north, densify_pts=21
    )
    left -= margin_m
    right += margin_m
    bottom -= margin_m
    top += margin_m
    width = max(1, int(math.ceil((right - left) / resolution_m)))
    height = max(1, int(math.ceil((top - bottom) / resolution_m)))
    dst_transform = from_bounds(left, bottom, right, top, width, height)
    return dst_crs, dst_transform, height, width


def _has_gcps(src) -> bool:
    gcps, gcrs = src.gcps
    return bool(gcps) and gcrs is not None


def _has_real_affine(src) -> bool:
    return src.crs is not None and not _is_identity_pixel_transform(src.transform)


def _warp_band_to_aoi_grid(
    path: str,
    dst_crs,
    dst_transform,
    dst_height: int,
    dst_width: int,
    resampling,
) -> np.ndarray:
    """Warp one band onto the AOI UTM grid as ``float32``.

    Three strategies, tried in order:

    1. **GCP path** (PC Sentinel-1 COGs): pass ``gcps=`` + ``src_crs=`` to ``reproject``
       (no ``src_transform``). This is the documented rasterio API for GCP files and
       avoids the ``GDALWarpOptions.Validate()`` crash on Linux GDAL.
    2. **Affine path** (CDSE / ASF .SAFE measurement TIFFs): normal ``src_transform`` +
       ``src_crs``.
    3. **osgeo.gdal.Warp fallback**: if both rasterio paths fail, use GDAL Python bindings
       directly — handles GCPs natively on every platform/version.
    """
    import logging
    import rasterio
    from rasterio.warp import reproject

    _log = logging.getLogger("pwtt.warp")
    fname = os.path.basename(path)
    out = np.empty((dst_height, dst_width), dtype=np.float32)

    # --- Strategy 1: GCP-based reproject (PC COGs) ---
    with rasterio.open(path) as src:
        if _has_gcps(src):
            gcps, gcrs = src.gcps
            arr = src.read(1).astype(np.float32, copy=False)
            try:
                reproject(
                    arr,
                    out,
                    gcps=gcps,
                    src_crs=gcrs,
                    dst_transform=dst_transform,
                    dst_crs=dst_crs,
                    resampling=resampling,
                )
                _log.info("%s: warped via GCP strategy", fname)
                return out
            except Exception as e1:
                _log.warning("%s: GCP strategy failed: %s", fname, e1)

    # --- Strategy 2: affine-based reproject (CDSE / ASF) ---
    with rasterio.open(path) as src:
        if _has_real_affine(src):
            arr = src.read(1).astype(np.float32, copy=False)
            try:
                reproject(
                    arr,
                    out,
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=dst_transform,
                    dst_crs=dst_crs,
                    resampling=resampling,
                )
                _log.info("%s: warped via affine strategy", fname)
                return out
            except Exception as e2:
                _log.warning("%s: affine strategy failed: %s", fname, e2)

    # --- Strategy 3: osgeo.gdal.Warp (universal fallback) ---
    _log.info("%s: falling back to osgeo.gdal.Warp", fname)
    return _gdal_warp_band_fallback(
        path, dst_crs, dst_transform, dst_height, dst_width, resampling
    )


def _gdal_warp_band_fallback(
    path: str,
    dst_crs,
    dst_transform,
    dst_height: int,
    dst_width: int,
    resampling,
) -> np.ndarray:
    """Last-resort warp using ``osgeo.gdal.Warp`` — works on every GDAL build."""
    from osgeo import gdal

    gdal.UseExceptions()

    _RESAMP_MAP = {
        0: gdal.GRA_NearestNeighbour,
        1: gdal.GRA_Bilinear,
        2: gdal.GRA_Cubic,
        3: gdal.GRA_CubicSpline,
        5: gdal.GRA_Average,
    }
    try:
        gra = _RESAMP_MAP.get(int(resampling), gdal.GRA_Bilinear)
    except Exception:
        gra = gdal.GRA_Bilinear

    dst_crs_wkt = dst_crs.to_wkt() if hasattr(dst_crs, "to_wkt") else str(dst_crs)
    gt = (
        dst_transform.c,
        dst_transform.a,
        dst_transform.b,
        dst_transform.f,
        dst_transform.d,
        dst_transform.e,
    )

    src_ds = gdal.Open(path, gdal.GA_ReadOnly)
    if src_ds is None:
        raise RuntimeError(f"GDAL cannot open {path}")

    mem_drv = gdal.GetDriverByName("MEM")
    dst_ds = mem_drv.Create("", dst_width, dst_height, 1, gdal.GDT_Float32)
    dst_ds.SetGeoTransform(gt)
    dst_ds.SetProjection(dst_crs_wkt)
    dst_ds.GetRasterBand(1).Fill(0)

    gdal.Warp(
        dst_ds,
        src_ds,
        resampleAlg=gra,
        dstSRS=dst_crs_wkt,
        width=dst_width,
        height=dst_height,
        outputBounds=(gt[0], gt[3] + gt[5] * dst_height, gt[0] + gt[1] * dst_width, gt[3]),
        outputBoundsSRS=dst_crs_wkt,
    )

    out = dst_ds.GetRasterBand(1).ReadAsArray().astype(np.float32, copy=False)
    dst_ds = None
    src_ds = None
    return out


def _short_progress_error(msg: str) -> str:
    if not msg:
        return "unknown error"
    one_line = " ".join(str(msg).split())
    if len(one_line) <= _PROGRESS_ERR_MAX_LEN:
        return one_line
    return one_line[: _PROGRESS_ERR_MAX_LEN - 1] + "…"


def _read_warp_vv_vh_pair(
    vv_path: str,
    vh_path: str,
    dst_crs,
    dst_transform,
    dst_height: int,
    dst_width: int,
    resampling,
    log,
) -> Optional[Tuple[np.ndarray, np.ndarray, dict, object, object]]:
    """Warp VV/VH to the AOI grid. Log errors and return ``None`` on failure."""
    try:
        vv = _warp_band_to_aoi_grid(
            vv_path, dst_crs, dst_transform, dst_height, dst_width, resampling
        )
        vh = _warp_band_to_aoi_grid(
            vh_path, dst_crs, dst_transform, dst_height, dst_width, resampling
        )
    except Exception as e:
        if log:
            log(f"warp VV/VH failed — {_short_progress_error(repr(e))}")
        return None
    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "width": dst_width,
        "height": dst_height,
        "count": 1,
        "crs": dst_crs,
        "transform": dst_transform,
    }
    return vv, vh, profile, dst_transform, dst_crs


def _add_months_dt(d: datetime, months: int) -> datetime:
    """Calendar add (same rule as openEO backend: clamp day to 28)."""
    m = d.month - 1 + months
    y = d.year + m // 12
    m = m % 12 + 1
    return datetime(y, m, min(d.day, 28))


def _circle_kernel(radius_m: float, pixel_size: float) -> np.ndarray:
    """Binary circle kernel in pixels."""
    r_px = int(np.ceil(radius_m / pixel_size))
    y, x = np.ogrid[-r_px : r_px + 1, -r_px : r_px + 1]
    return ((x * x + y * y) <= (radius_m / pixel_size) ** 2).astype(np.float64)


def _settings_local_source() -> str:
    s = QgsSettings()
    s.beginGroup("PWTT")
    src = (s.value("local_data_source", LOCAL_SOURCE_CDSE) or LOCAL_SOURCE_CDSE).strip().lower()
    s.endGroup()
    if src not in (LOCAL_SOURCE_CDSE, LOCAL_SOURCE_ASF, LOCAL_SOURCE_PC):
        return LOCAL_SOURCE_CDSE
    return src


class LocalBackend(PWTTBackend):
    @property
    def name(self):
        return "Local Processing (Experimental)"

    @property
    def id(self):
        return "local"

    def check_dependencies(self):
        from . import deps

        src = _settings_local_source()
        missing, pip = deps.local_backend_missing(src)
        if missing:
            if pip:
                return False, f"Local backend ({src}) requires: pip install {' '.join(pip)}"
            return False, f"Local backend missing: {', '.join(missing)}"
        return True, ""

    def authenticate(self, credentials: dict) -> bool:
        source = (credentials.get("source") or LOCAL_SOURCE_CDSE).strip().lower()
        if source not in (LOCAL_SOURCE_CDSE, LOCAL_SOURCE_ASF, LOCAL_SOURCE_PC):
            source = LOCAL_SOURCE_CDSE
        self._data_source = source

        if source == LOCAL_SOURCE_CDSE:
            user = (credentials.get("username") or "").strip()
            password = credentials.get("password") or ""
            if not user:
                raise ValueError("CDSE username is required.")
            if not password:
                raise ValueError("CDSE password is required.")
            try:
                self._token = get_token(user, password)
            except Exception as e:
                raise RuntimeError(f"CDSE authentication failed: {e}") from e
            self._asf_session = None
            self._pc_client = None
            return True

        if source == LOCAL_SOURCE_ASF:
            from .asf_downloader import authenticate_asf

            user = (credentials.get("earthdata_username") or credentials.get("username") or "").strip()
            password = credentials.get("earthdata_password") or credentials.get("password") or ""
            if not user:
                raise ValueError("Earthdata username is required for ASF.")
            if not password:
                raise ValueError("Earthdata password is required for ASF.")
            try:
                self._asf_session = authenticate_asf(user, password)
            except Exception as e:
                raise RuntimeError(
                    "ASF / Earthdata authentication failed. "
                    "Check Earthdata credentials and account access for ASF downloads. "
                    f"Details: {e}"
                ) from e
            self._token = None
            self._pc_client = None
            return True

        if source == LOCAL_SOURCE_PC:
            from .pc_downloader import authenticate_pc

            key = (credentials.get("pc_subscription_key") or "").strip() or None
            try:
                self._pc_client = authenticate_pc(key)
            except Exception as e:
                raise RuntimeError(
                    "Planetary Computer client initialization failed. "
                    "Check network access and subscription key (if provided). "
                    f"Details: {e}"
                ) from e
            self._token = None
            self._asf_session = None
            return True

        return False

    def run(
        self,
        aoi_wkt: str,
        war_start: str,
        inference_start: str,
        pre_interval: int,
        post_interval: int,
        output_path: str,
        progress_callback=None,
        include_footprints: bool = False,
        footprints_path: Optional[str] = None,
        remote_job_id: Optional[str] = None,
        damage_threshold: float = 3.3,
        gee_viz: bool = False,
        save_timeseries: bool = True,
        job_id: Optional[str] = None,
    ) -> str:
        source = getattr(self, "_data_source", LOCAL_SOURCE_CDSE)
        if source == LOCAL_SOURCE_CDSE and not getattr(self, "_token", None):
            raise ValueError("Not authenticated. Call authenticate() first.")
        if source == LOCAL_SOURCE_ASF and not getattr(self, "_asf_session", None):
            raise ValueError("Not authenticated. Call authenticate() first.")
        if source == LOCAL_SOURCE_PC and not getattr(self, "_pc_client", None):
            raise ValueError("Not authenticated. Call authenticate() first.")

        import rasterio
        from rasterio.warp import Resampling

        bbox = wkt_to_bbox(aoi_wkt)
        if not bbox:
            raise ValueError("Invalid AOI WKT")
        west, south, east, north = bbox
        dst_crs, dst_transform, dst_height, dst_width = _aoi_utm_grid(
            west, south, east, north
        )
        war_d = datetime.strptime(war_start[:10], "%Y-%m-%d")
        inf_d = datetime.strptime(inference_start[:10], "%Y-%m-%d")
        pre_start = _add_months_dt(war_d, -pre_interval).strftime("%Y-%m-%d")
        post_end = _add_months_dt(inf_d, post_interval).strftime("%Y-%m-%d")

        cache_dir = os.path.join(os.path.dirname(output_path), ".pwtt_cache")
        os.makedirs(cache_dir, exist_ok=True)

        # QgsTask progress must not be reset to 0 on every log line (search/download are noisy).
        _progress = [0]

        def emit(pct: int, msg: str = ""):
            pct = max(0, min(100, int(pct)))
            _progress[0] = max(_progress[0], pct)
            if progress_callback:
                progress_callback(_progress[0], msg)

        def job_log(msg: str):
            if progress_callback:
                progress_callback(_progress[0], msg)

        job_log(
            f"Local processing: source={source}, cache_dir={cache_dir}, "
            f"output_tif={output_path}"
        )
        job_log(
            f"Local: AOI WGS84 bbox west={west:.5f} south={south:.5f} "
            f"east={east:.5f} north={north:.5f}"
        )
        job_log(
            f"Local: AOI warp ~{_AOI_RESOLUTION_M} m UTM, {dst_width}×{dst_height} px, "
            f"CRS {dst_crs} (per-scene warp; pre/post loaded sequentially)"
        )
        job_log(
            f"Local: pre window {pre_start} … {war_start[:10]}, "
            f"post window {inference_start[:10]} … {post_end}"
        )

        max_per_period = _max_scenes_per_period()
        self.run_metadata = {
            "collection": "SENTINEL-1 IW_GRDH_1S",
            "data_source": source,
            "sar_backscatter": "sigma0-linear-grd (openEO-aligned; no Lee/log)",
            "pre_period": {"start": pre_start, "end": war_start},
            "post_period": {"start": inference_start, "end": post_end},
            "bbox": [west, south, east, north],
            "pre_scenes_found": [],
            "post_scenes_found": [],
            "pre_scenes_used": [],
            "post_scenes_used": [],
            "offline_scenes": [],
            "local_max_scenes_per_period": max_per_period,
        }

        emit(5, "Searching pre-war/event products…")

        if source == LOCAL_SOURCE_CDSE:
            pre_products = search_s1_grd(
                self._token,
                aoi_wkt,
                pre_start,
                war_start,
                max_results=80,
                log=job_log,
            )
            emit(10, "Searching post-war/event products…")
            post_products = search_s1_grd(
                self._token,
                aoi_wkt,
                inference_start,
                post_end,
                max_results=80,
                log=job_log,
            )
        elif source == LOCAL_SOURCE_ASF:
            from .asf_downloader import search_s1_grd_asf

            emit(5, "Searching pre-war/event products…")
            pre_products = search_s1_grd_asf(
                self._asf_session,
                aoi_wkt,
                pre_start,
                war_start,
                max_results=80,
                log=job_log,
            )
            emit(10, "Searching post-war/event products…")
            post_products = search_s1_grd_asf(
                self._asf_session,
                aoi_wkt,
                inference_start,
                post_end,
                max_results=80,
                log=job_log,
            )
        else:
            from .pc_downloader import search_s1_grd_pc

            emit(5, "Searching pre-war/event products…")
            pre_products = search_s1_grd_pc(
                self._pc_client,
                aoi_wkt,
                pre_start,
                war_start,
                max_results=80,
                log=job_log,
            )
            emit(10, "Searching post-war/event products…")
            post_products = search_s1_grd_pc(
                self._pc_client,
                aoi_wkt,
                inference_start,
                post_end,
                max_results=80,
                log=job_log,
            )

        if not pre_products or not post_products:
            if source == LOCAL_SOURCE_CDSE:
                raise RuntimeError(
                    "No Sentinel-1 GRD products found on CDSE for this AOI/date range. "
                    "Try a wider date range or switch local source to ASF/Planetary Computer."
                )
            if source == LOCAL_SOURCE_ASF:
                raise RuntimeError(
                    "No Sentinel-1 GRD products found on ASF for this AOI/date range. "
                    "Confirm Earthdata access and try a wider date range."
                )
            raise RuntimeError(
                "No Sentinel-1 GRD products found on Planetary Computer for this AOI/date range. "
                "Try a wider date range or use CDSE/ASF source."
            )

        def _scene_summary(prod):
            cd = prod.get("ContentDate") or {}
            name = prod.get("Name", "")
            return {
                "id": prod.get("Id", ""),
                "name": name,
                "date": (cd.get("Start") or "")[:19],
                "online": prod.get("Online", True),
                "orbit": _parse_relative_orbit_from_name(name),
            }

        self.run_metadata["pre_scenes_found"] = [_scene_summary(p) for p in pre_products]
        self.run_metadata["post_scenes_found"] = [_scene_summary(p) for p in post_products]

        def _preview_names(products, n=6):
            parts = []
            for p in products[:n]:
                parts.append(str(p.get("Name") or p.get("Id") or "?"))
            tail = f" (+{len(products) - n} more)" if len(products) > n else ""
            return ", ".join(parts) + tail

        job_log(f"Local: pre search — {_preview_names(pre_products)}")
        job_log(f"Local: post search — {_preview_names(post_products)}")
        job_log(
            f"Local: will try up to {max_per_period} pre + {max_per_period} post scenes "
            f"(openEO-aligned σ⁰ linear, pooled t; QgsSettings PWTT/local_max_scenes_per_period)"
        )
        pre_arrays = []
        post_arrays = []
        triggered_orders = [0]  # mutable container for thread-safe increment
        offline_product_ids = []
        download_failures = []
        _lock = threading.Lock()

        def load_products_safe(products, arrays_list, label, used_key):
            loaded = 0
            for prod in products:
                if loaded >= max_per_period:
                    break
                emit(
                    min(54, 12 + int(42 * (loaded + 1) / max_per_period)),
                    f"{label} {loaded + 1}/{max_per_period}…",
                )
                pid, name = prod["Id"], prod["Name"]
                cd = prod.get("ContentDate") or {}
                scene_date = (cd.get("Start") or "")[:19]

                if source == LOCAL_SOURCE_CDSE:
                    safe_dir = download_product(
                        self._token,
                        pid,
                        name,
                        cache_dir,
                        wait_for_offline=False,
                        log=job_log,
                    )
                    if safe_dir is None:
                        with _lock:
                            triggered_orders[0] += 1
                            offline_product_ids.append(pid)
                            self.run_metadata["offline_scenes"].append(
                                {"id": pid, "name": name, "date": scene_date}
                            )
                        job_log(f"{label}: {name} is offline, staging order triggered…")
                        continue
                else:
                    from .asf_downloader import download_product_asf

                    err_msg = None
                    try:
                        safe_dir = download_product_asf(
                            self._asf_session, prod, cache_dir, log=job_log
                        )
                    except Exception as dl_err:
                        safe_dir = None
                        err_msg = str(dl_err)
                        with _lock:
                            download_failures.append({"name": name, "error": err_msg})
                    if not safe_dir:
                        if err_msg is None:
                            err_msg = "no SAFE/ZIP in cache after download (unexpected)"
                            with _lock:
                                download_failures.append({"name": name, "error": err_msg})
                        detail = _short_progress_error(err_msg)
                        job_log(f"{label}: skip {name} — {detail}")
                        continue

                vv_path, vh_path = find_vv_vh_in_safe(safe_dir)
                if not vv_path or not vh_path:
                    job_log(
                        f"{label}: no VV/VH tifs under measurement/ in {safe_dir} "
                        f"(skip «{name}»)"
                    )
                    continue
                job_log(
                    f"{label}: VV/VH rasters — {vv_path} | {vh_path}"
                )
                result = _read_warp_vv_vh_pair(
                    vv_path,
                    vh_path,
                    dst_crs,
                    dst_transform,
                    dst_height,
                    dst_width,
                    Resampling.bilinear,
                    job_log,
                )
                if result is None:
                    job_log(f"{label}: skip {name} (AOI warp failed)…")
                    continue
                if source == LOCAL_SOURCE_CDSE:
                    from .downloader import remove_product_zip

                    remove_product_zip(name, cache_dir, log=job_log)
                elif source == LOCAL_SOURCE_ASF:
                    from .asf_downloader import remove_zips_for_extracted_safe

                    remove_zips_for_extracted_safe(safe_dir, log=job_log)
                arrays_list.append(result)
                vv_arr, vh_arr = result[0], result[1]
                try:
                    vv_mean = float(np.nanmean(vv_arr))
                except Exception:
                    vv_mean = float("nan")
                try:
                    vh_mean = float(np.nanmean(vh_arr))
                except Exception:
                    vh_mean = float("nan")
                with _lock:
                    self.run_metadata[used_key].append({
                        "id": pid,
                        "name": name,
                        "date": scene_date,
                        "orbit": _parse_relative_orbit_from_name(name),
                        "VV_mean": vv_mean,
                        "VH_mean": vh_mean,
                    })
                loaded += 1
                job_log(
                    f"{label}: accepted scene «{name}» ({scene_date}) "
                    f"[{loaded}/{max_per_period}]"
                )

        def load_products_pc(products, arrays_list, label, used_key):
            from .pc_downloader import download_pc_vv_vh

            loaded = 0
            for prod in products:
                if loaded >= max_per_period:
                    break
                emit(
                    min(54, 12 + int(42 * (loaded + 1) / max_per_period)),
                    f"{label} {loaded + 1}/{max_per_period}…",
                )
                pid, name = prod["Id"], prod["Name"]
                cd = prod.get("ContentDate") or {}
                scene_date = (cd.get("Start") or "")[:19]
                subdir = os.path.join(cache_dir, "pc_" + "".join(c if c.isalnum() else "_" for c in pid)[:120])
                dl_err = None
                try:
                    vv_path, vh_path = download_pc_vv_vh(prod, subdir, log=job_log)
                except Exception as e:
                    dl_err = e
                    vv_path, vh_path = None, None
                    with _lock:
                        download_failures.append({"name": name, "error": str(dl_err)})
                if not vv_path or not vh_path:
                    if dl_err is not None:
                        detail = _short_progress_error(str(dl_err))
                        job_log(f"{label}: skip {name} — {detail}")
                    else:
                        job_log(f"{label}: skip {name} (no VV/VH assets)…")
                    continue
                job_log(f"{label}: PC COGs ready — {vv_path} | {vh_path}")
                result = _read_warp_vv_vh_pair(
                    vv_path,
                    vh_path,
                    dst_crs,
                    dst_transform,
                    dst_height,
                    dst_width,
                    Resampling.bilinear,
                    job_log,
                )
                if result is None:
                    job_log(f"{label}: skip {name} (AOI warp failed)…")
                    continue
                arrays_list.append(result)
                vv_arr, vh_arr = result[0], result[1]
                try:
                    vv_mean = float(np.nanmean(vv_arr))
                except Exception:
                    vv_mean = float("nan")
                try:
                    vh_mean = float(np.nanmean(vh_arr))
                except Exception:
                    vh_mean = float("nan")
                with _lock:
                    self.run_metadata[used_key].append({
                        "id": pid,
                        "name": name,
                        "date": scene_date,
                        "orbit": _parse_relative_orbit_from_name(name),
                        "VV_mean": vv_mean,
                        "VH_mean": vh_mean,
                    })
                loaded += 1
                job_log(
                    f"{label}: accepted PC scene «{name}» ({scene_date}) "
                    f"[{loaded}/{max_per_period}]"
                )

        emit(12, "Downloading pre- and post-war/event scenes…")

        loader = load_products_pc if source == LOCAL_SOURCE_PC else load_products_safe

        # Serial pre/post: each scene warp peaks at ~1–2 full GRD bands + small AOI;
        # parallel load previously paired with MemoryError / GDAL thread issues in QGIS.
        loader(pre_products, pre_arrays, "Pre", "pre_scenes_used")
        loader(post_products, post_arrays, "Post", "post_scenes_used")

        job_log(
            f"Local: downloads finished — pre stacks={len(pre_arrays)}, "
            f"post stacks={len(post_arrays)}"
        )

        if not pre_arrays or not post_arrays:
            if source == LOCAL_SOURCE_CDSE and triggered_orders[0] > 0:
                from .base_backend import ProductsOfflineError

                raise ProductsOfflineError(
                    f"All available products are in cold storage. "
                    f"Staging orders have been triggered for {triggered_orders[0]} product(s). "
                    f"Will auto-check and resume when products become available.",
                    product_ids=offline_product_ids,
                    offline_scenes=list(self.run_metadata.get("offline_scenes", [])),
                )
            if source == LOCAL_SOURCE_ASF:
                detail = ""
                if download_failures:
                    detail = " Errors: " + "; ".join(
                        f"{f['name']}: {f['error']}" for f in download_failures[:3]
                    )
                raise RuntimeError(
                    "Could not load VV/VH data from ASF products. "
                    "Try a different date range or verify Earthdata/ASF access."
                    + detail
                )
            if source == LOCAL_SOURCE_PC:
                detail = ""
                if download_failures:
                    detail = " Errors: " + "; ".join(
                        f"{f['name']}: {f['error']}" for f in download_failures[:3]
                    )
                raise RuntimeError(
                    "Could not load VV/VH assets from Planetary Computer STAC items. "
                    "Try a different date range or switch source."
                    + detail
                )
            raise RuntimeError(
                "Could not load VV/VH data from any CDSE product. "
                "All products may be offline; try again later or switch source."
            )

        ref_vv, ref_vh, ref_profile, ref_transform, ref_crs = pre_arrays[0]
        height, width = ref_vv.shape
        pixel_size = abs(ref_transform.a)
        job_log(
            f"Local: reference grid {width}×{height} px, "
            f"~{pixel_size:.2f} m, CRS {ref_crs}"
        )

        def _accumulate_band(arrays, band_idx: int):
            """band_idx 0 = VV, 1 = VH. σ⁰ linear on warped GRD (matches openEO sar_backscatter)."""
            mean, m2, n = welford_init((height, width))
            for tpl in arrays:
                vv, vh = tpl[0], tpl[1]
                x = vv if band_idx == 0 else vh
                mean, m2, n = welford_update(mean, m2, n, x)
            return mean, welford_sample_variance(m2, n), n.astype(np.float64)

        job_log(
            f"Local: σ⁰ linear (no Lee/log), Welford mean/var — "
            f"pre {len(pre_arrays)} scenes, post {len(post_arrays)} scenes"
        )
        emit(55, "Computing t-test…")

        eps = 1e-12
        pre_mean_vv, pre_var_vv, pre_n_vv = _accumulate_band(pre_arrays, 0)
        pre_mean_vh, pre_var_vh, pre_n_vh = _accumulate_band(pre_arrays, 1)
        post_mean_vv, post_var_vv, post_n_vv = _accumulate_band(post_arrays, 0)
        post_mean_vh, post_var_vh, post_n_vh = _accumulate_band(post_arrays, 1)

        df_vv = np.maximum(pre_n_vv + post_n_vv - 2.0, 1.0)
        df_vh = np.maximum(pre_n_vh + post_n_vh - 2.0, 1.0)
        pooled_var_vv = (
            pre_var_vv * np.maximum(pre_n_vv - 1.0, 0.0)
            + post_var_vv * np.maximum(post_n_vv - 1.0, 0.0)
        ) / df_vv
        pooled_var_vh = (
            pre_var_vh * np.maximum(pre_n_vh - 1.0, 0.0)
            + post_var_vh * np.maximum(post_n_vh - 1.0, 0.0)
        ) / df_vh
        pooled_sd_vv = np.sqrt(np.maximum(pooled_var_vv, 0.0)) + eps
        pooled_sd_vh = np.sqrt(np.maximum(pooled_var_vh, 0.0)) + eps
        inv_sqrt_n_vv = np.sqrt(
            1.0 / np.maximum(pre_n_vv, 1.0) + 1.0 / np.maximum(post_n_vv, 1.0)
        )
        inv_sqrt_n_vh = np.sqrt(
            1.0 / np.maximum(pre_n_vh, 1.0) + 1.0 / np.maximum(post_n_vh, 1.0)
        )
        se_vv = pooled_sd_vv * inv_sqrt_n_vv + eps
        se_vh = pooled_sd_vh * inv_sqrt_n_vh + eps
        t_vv = np.abs(post_mean_vv - pre_mean_vv) / se_vv
        t_vh = np.abs(post_mean_vh - pre_mean_vh) / se_vh
        max_change = np.maximum(t_vv, t_vh)
        job_log(
            f"Local: pooled t (openEO-style) — scenes pre={len(pre_arrays)} post={len(post_arrays)}, "
            f"damage threshold={damage_threshold}"
        )

        p_value = openeo_style_p_value_bound(max_change)

        emit(70, "Post-processing…")
        job_log("Local: 50/100/150 m circular means on max |t| (no extra Gaussian; openEO-aligned)")

        def _mean_kernel(radius_m):
            k = _circle_kernel(radius_m, pixel_size)
            k = k / (k.sum() + 1e-12)
            return k

        k50 = convolve2d_edge(max_change, _mean_kernel(50))
        k100 = convolve2d_edge(max_change, _mean_kernel(100))
        k150 = convolve2d_edge(max_change, _mean_kernel(150))
        t_statistic = (max_change + k50 + k100 + k150) / 4.0
        damage = (t_statistic > float(damage_threshold)).astype(np.float32)

        emit(90, "Writing GeoTIFF…")
        job_log(
            f"Local: write GeoTIFF bands (1=t, 2=damage mask, 3=p-value) LZW → {output_path}"
        )
        out_profile = ref_profile.copy()
        out_profile.update(
            dtype=rasterio.float32,
            count=3,
            compress="lzw",
            nodata=-9999,
        )
        with rasterio.open(output_path, "w", **out_profile) as dst:
            dst.write(t_statistic.astype(np.float32), 1)
            dst.write(damage, 2)
            dst.write(p_value.astype(np.float32), 3)

        self.run_metadata["output_size_bytes"] = os.path.getsize(output_path)
        self.run_metadata["output_crs"] = str(ref_crs)
        self.run_metadata["output_pixel_size_m"] = round(pixel_size, 2)
        self.run_metadata["output_shape"] = [height, width]
        self.run_metadata["pre_scenes_count"] = len(pre_arrays)
        self.run_metadata["post_scenes_count"] = len(post_arrays)

        if save_timeseries:
            try:
                self._write_local_timeseries_sidecar(
                    output_path=output_path,
                    aoi_wkt=aoi_wkt,
                    war_start=war_start,
                    inference_start=inference_start,
                    pre_interval=pre_interval,
                    post_interval=post_interval,
                    job_id=job_id or "",
                )
            except Exception as ts_err:
                job_log(f"Local: time-series sidecar skipped — {ts_err}")

        emit(95, "Done.")
        return output_path

    def _write_local_timeseries_sidecar(
        self,
        output_path: str,
        aoi_wkt: str,
        war_start: str,
        inference_start: str,
        pre_interval: int,
        post_interval: int,
        job_id: str,
    ) -> None:
        """Pool pre-war/event means, compute VV/VH z-scores per scene, write sidecars.

        Local backend pools all orbits into one composite, so normalization here
        is pooled (not per-orbit) — the JSON records this explicitly.
        """
        pre_used = list(self.run_metadata.get("pre_scenes_used") or [])
        post_used = list(self.run_metadata.get("post_scenes_used") or [])
        if not pre_used and not post_used:
            return

        def _finite(xs):
            out = []
            for x in xs:
                try:
                    if x is None:
                        continue
                    fx = float(x)
                    if fx == fx and fx not in (float("inf"), float("-inf")):
                        out.append(fx)
                except Exception:
                    continue
            return out

        pre_vv = _finite(s.get("VV_mean") for s in pre_used)
        pre_vh = _finite(s.get("VH_mean") for s in pre_used)
        if len(pre_vv) < 2 or len(pre_vh) < 2:
            return

        pre_vv_mean = float(np.mean(pre_vv))
        pre_vh_mean = float(np.mean(pre_vh))
        pre_vv_std = float(np.std(pre_vv, ddof=1))
        pre_vh_std = float(np.std(pre_vh, ddof=1))

        def _z(x, mu, sd):
            if x is None or sd is None or sd == 0.0 or sd != sd:
                return None
            try:
                fx = float(x)
            except Exception:
                return None
            if fx != fx:
                return None
            return (fx - mu) / sd

        def _entry(scene, period):
            vv_z = _z(scene.get("VV_mean"), pre_vv_mean, pre_vv_std)
            vh_z = _z(scene.get("VH_mean"), pre_vh_mean, pre_vh_std)
            return {
                "date": scene.get("date") or "",
                "orbit": scene.get("orbit"),
                "pass": None,
                "VV_z": vv_z,
                "VH_z": vh_z,
                "period": period,
            }

        series = []
        for s in pre_used:
            series.append(_entry(s, "pre"))
        for s in post_used:
            series.append(_entry(s, "post"))
        series.sort(key=lambda e: (e.get("date") or ""))

        from . import timeseries_sidecar
        payload = timeseries_sidecar.build_sidecar(
            job_id=job_id,
            backend=self.id,
            aoi_wkt=aoi_wkt,
            war_start=war_start,
            inference_start=inference_start,
            pre_interval_months=pre_interval,
            post_interval_months=post_interval,
            normalization=(
                "pooled z-score vs pre-war/event baseline (Local backend pools orbits; "
                "mean/std over σ⁰ AOI means)"
            ),
            series=series,
        )
        timeseries_sidecar.write_sidecars(output_path, payload)
