"""Per-tile LiDAR data quality: ground point density and acquisition date.

One JSON sidecar per tile (output/quality/{basename}.json) is written by the
pipeline, copied into the inventory (index_tiles.json, `quality` table) and
kept on the lightweight machines: the quality inset of the PDF export reads it
from the local disk, whether the upstream is reachable or not.

numpy is only imported by the computation: the lightweight image (without
numpy) reads and writes the sidecars.
"""

import json
import logging
import os
import tempfile
from datetime import datetime, timedelta, timezone
from pathlib import Path

logger = logging.getLogger("lidar")

QUALITY_VERSION = 1
QUALITY_DIRNAME = "quality"
DENSITY_CELL_M = 50.0      # density grid cell size (20 × 20 per tile)
EMPTY_PIXEL_M = 0.2        # DTM pixel: share without ground point ≈ interpolated share

_GPS_EPOCH = datetime(1980, 1, 6, tzinfo=timezone.utc)
_GPS_ADJUSTED_OFFSET = 1e9  # standard adjusted GPS time = GPS seconds − 1e9


def gps_adjusted_to_date(t):
    """Standard adjusted GPS time → ISO date (UTC; leap seconds ignored:
    accurate to the day)."""
    return (_GPS_EPOCH + timedelta(seconds=float(t) + _GPS_ADJUSTED_OFFSET)).date().isoformat()


def compute_quality(x, y, gps_time, bounds, gps_adjusted=True, header_date=None):
    """Measure the quality of a tile from its ground points.

    Args:
        x, y: L93 coordinates of the ground points.
        gps_time: GPS time of the points (or None).
        bounds: nominal footprint (min_x, min_y, max_x, max_y); points
            outside it (edge buffer) are ignored.
        gps_adjusted: True if gps_time is standard adjusted GPS time.
        header_date: ISO date from the LAS header (date fallback).

    Returns:
        JSON-serializable dict: version, ground_density (pts/m²),
        density_grid (pts/m² per 50 m cell, rows north to south),
        empty_fraction (share of 0.2 m pixels without a ground point),
        acq_start / acq_end (ISO dates) and acq_source ("gps", "header" or None).
    """
    import numpy as np

    min_x, min_y, max_x, max_y = bounds
    x = np.asarray(x, dtype=np.float64)
    y = np.asarray(y, dtype=np.float64)
    inside = (x >= min_x) & (x < max_x) & (y >= min_y) & (y < max_y)
    t = None
    if gps_time is not None and gps_adjusted:
        t = np.asarray(gps_time, dtype=np.float64)
        t = t[inside] if len(t) == len(x) else None
    x, y = x[inside], y[inside]
    n = len(x)
    width, height = max_x - min_x, max_y - min_y

    nx = int(round(width / DENSITY_CELL_M))
    ny = int(round(height / DENSITY_CELL_M))
    cx = np.clip(((x - min_x) / DENSITY_CELL_M).astype(np.int64), 0, nx - 1)
    cy = np.clip(((max_y - y) / DENSITY_CELL_M).astype(np.int64), 0, ny - 1)
    counts = np.bincount(cy * nx + cx, minlength=nx * ny).reshape(ny, nx)
    grid = np.round(counts / (DENSITY_CELL_M ** 2), 2)

    px = int(round(width / EMPTY_PIXEL_M))
    py = int(round(height / EMPTY_PIXEL_M))
    occupied = np.zeros(px * py, dtype=bool)
    if n:
        ix = np.clip(((x - min_x) / EMPTY_PIXEL_M).astype(np.int64), 0, px - 1)
        iy = np.clip(((max_y - y) / EMPTY_PIXEL_M).astype(np.int64), 0, py - 1)
        occupied[iy * px + ix] = True
    empty_fraction = 1.0 - float(occupied.sum()) / (px * py)

    if t is not None and len(t):
        start, end = gps_adjusted_to_date(t.min()), gps_adjusted_to_date(t.max())
        source = "gps"
    elif header_date:
        start = end = header_date
        source = "header"
    else:
        start = end = source = None

    return {
        "version": QUALITY_VERSION,
        "ground_density": round(n / (width * height), 6),
        "density_grid": grid.tolist(),
        "empty_fraction": round(empty_fraction, 4),
        "acq_start": start,
        "acq_end": end,
        "acq_source": source,
    }


def quality_path(output_dir, basename):
    """Path of a tile's quality sidecar."""
    return Path(output_dir) / QUALITY_DIRNAME / f"{basename}.json"


def write_quality(output_dir, basename, data):
    """Write the sidecar (atomically); False if the content is already identical."""
    path = quality_path(output_dir, basename)
    payload = json.dumps(data, ensure_ascii=False, sort_keys=True)
    try:
        if path.read_text(encoding="utf-8") == payload:
            return False
    except OSError:
        pass
    path.parent.mkdir(parents=True, exist_ok=True)
    fd, tmp = tempfile.mkstemp(dir=str(path.parent), suffix=".tmp")
    try:
        with os.fdopen(fd, "w", encoding="utf-8") as f:
            f.write(payload)
        os.replace(tmp, path)
    except OSError:
        try:
            os.unlink(tmp)
        except OSError:
            pass
        raise
    return True


def read_quality(output_dir, basename):
    """A tile's sidecar, or None (missing, unreadable, other version)."""
    try:
        data = json.loads(quality_path(output_dir, basename).read_text(encoding="utf-8"))
    except (OSError, ValueError):
        return None
    if not isinstance(data, dict) or data.get("version") != QUALITY_VERSION:
        return None
    return data


def load_quality_table(output_dir):
    """All valid sidecars: {basename: sidecar}."""
    folder = Path(output_dir) / QUALITY_DIRNAME
    table = {}
    if not folder.is_dir():
        return table
    for f in sorted(folder.glob("*.json")):
        data = read_quality(output_dir, f.stem)
        if data is not None:
            table[f.stem] = data
    return table


_LAZ_EXTS = (".copc.laz", ".copc.las", ".laz", ".las")


def _laz_basename(path):
    """Tile name without its LiDAR extension (.copc.laz included)."""
    name = Path(path).name
    for ext in _LAZ_EXTS:
        if name.lower().endswith(ext):
            return name[:-len(ext)]
    return Path(path).stem


def _cell_bounds_of(basename):
    """Nominal 1 km footprint of an LHD tile, or None outside the LHD pattern."""
    from .index import parse_basename_coords
    coords = parse_basename_coords(basename)
    if coords is None:
        return None
    col, row = coords
    return (col * 1000.0, (row - 1) * 1000.0, (col + 1) * 1000.0, row * 1000.0)


def quality_from_las(las_path, bounds, codes=None):
    """Quality measured from a LAS/LAZ; None on failure (never raises).

    codes: classes to keep (returns ≥ 1); None = file already filtered to ground.
    """
    try:
        import laspy
        import numpy as np
        las = laspy.read(str(las_path))
        keep = np.ones(len(las.points), dtype=bool)
        if codes is not None:
            keep = ((np.asarray(las.return_number) >= 1)
                    & (np.asarray(las.number_of_returns) >= 1)
                    & np.isin(np.asarray(las.classification), np.asarray(sorted(codes))))
        try:
            gps = np.asarray(las.gps_time, dtype=np.float64)[keep]
        except AttributeError:
            gps = None
        adjusted = las.header.global_encoding.gps_time_type == laspy.header.GpsTimeType.STANDARD
        created = las.header.creation_date
        return compute_quality(np.asarray(las.x)[keep], np.asarray(las.y)[keep], gps,
                               bounds, gps_adjusted=adjusted,
                               header_date=created.isoformat() if created else None)
    except Exception as e:  # noqa: BLE001 — best-effort quality
        logger.warning(f"  Could not measure quality ({Path(las_path).name}): {e}")
        return None


def ensure_quality(las_path, basename, output_dir, codes=None):
    """Ensure a tile's quality sidecar exists (computed if missing or outdated).

    Returns:
        True if a valid sidecar exists on return.
    """
    if read_quality(output_dir, basename) is not None:
        return True
    bounds = _cell_bounds_of(basename)
    if bounds is None:
        return False
    data = quality_from_las(las_path, bounds, codes)
    if data is None:
        return False
    try:
        write_quality(output_dir, basename, data)
    except OSError as e:
        logger.warning(f"  Could not write the quality sidecar: {e}")
        return False
    return True


def backfill_quality(input_dir, output_dir, codes=(2,)):
    """Backfill: quality sidecar for every LHD LAZ in input_dir that lacks one.

    Returns:
        Number of sidecars written.
    """
    files = sorted(p for p in Path(input_dir).iterdir()
                   if p.is_file() and p.name.lower().endswith(_LAZ_EXTS))
    written = 0
    for i, laz in enumerate(files, 1):
        base = _laz_basename(laz)
        if read_quality(output_dir, base) is not None or _cell_bounds_of(base) is None:
            continue
        logger.info(f"Quality [{i}/{len(files)}] {base}")
        if ensure_quality(laz, base, output_dir, codes):
            written += 1
    return written