"""DTM generation from classified LiDAR point clouds.

Handles ground classification via PDAL (IGN supplier pre-classification,
SMRF or CSF) and DTM rasterisation
using scipy binned_statistic_2d. Gaps without LiDAR data (common in
complex/rocky terrain) are filled with a terrain-aware interpolation so the
DTM stays continuous.
"""

import json
import logging
import subprocess
from pathlib import Path

import numpy as np
import rasterio
from rasterio.transform import from_bounds
from scipy.stats import binned_statistic_2d

logger = logging.getLogger("lidar")

# Classes LAS exploitables de la pré-classification LiDAR HD (noms → codes)
IGN_CLASS_NAMES = {
    "sol": 2,
    "unclassified": 1,
    "non-classe": 1,
    "eau": 9,
    "virtuel": 66,
    "pont": 17,
    "sursol": 64,
}


# Calage vertical des faisceaux de vol (strip alignment). Une tuile LiDAR HD
# est couverte par plusieurs passes d'acquisition, chacune portée par un ou
# deux PointSourceId (faisceaux). Les passes sont bien alignées horizontalement
# mais certains faisceaux portent un biais vertical de quelques cm (mesuré
# jusqu'à ~5 cm sur LHD_FXX_1000_6882 : les deux faisceaux d'une passe écartés
# de ±2,5 cm). À 0,2 m/px ces écarts créent des marches et du bruit aux
# coutures des zones de recouvrement. On mesure l'offset robuste de chaque
# faisceau sur les points sol de la tuile elle-même (les PointSourceId changent
# selon la campagne, rien n'est codé en dur) et on le retranche avant
# rastérisation. Les offsets dérivent le long d'une ligne de vol (signe inversé
# entre tuiles voisines mesuré) : le calcul est donc par tuile, jamais global.
#
# 2ᵉ passe — gigue intra-faisceau : le décalage vertical peut aussi varier au
# fil d'une MÊME passe (vibration capteur / bruit haute fréquence de la
# trajectoire) : les lignes de balayage successives d'un faisceau apparaissent
# alors décalées verticalement de façon aléatoire. Chaque faisceau est découpé
# en fenêtres de temps GPS, l'offset robuste de chaque fenêtre est mesuré contre
# la surface médiane des autres faisceaux (même maille 1 m), la série est lissée
# (médiane glissante) puis interpolée au temps GPS de chaque point. Requiert la
# dimension gps_time (ignorée silencieusement sinon).
STRIP_ALIGN_VERSION = 2
STRIP_ALIGN_THRESHOLD = 0.005  # m : écart mini pour corriger un faisceau (0,5 cm)
STRIP_ALIGN_CELL = 1.0         # m : maille de comparaison des faisceaux
STRIP_ALIGN_MIN_SHARED = 500   # cellules sol communes mini pour valider un offset
STRIP_JITTER_BIN = 0.1         # s : durée d'une fenêtre de temps GPS (gigue)
STRIP_JITTER_SMOOTH = 5        # fenêtres : largeur de la médiane glissante
STRIP_JITTER_MIN_CELLS = 40    # cellules sol partagées mini pour valider une fenêtre
STRIP_JITTER_MAX = 0.10        # m : amplitude maxi d'une correction de gigue (garde-fou)

# Mémo des offsets par fichier : la classification est partagée entre
# résolutions, le même LAS sol est rasterisé à 0,5 m puis 0,2 m.
_STRIP_OFFSETS_CACHE = {}
_STRIP_JITTER_CACHE = {}


def _strip_surface_grid(x, y, z, inv, n_sources, cell=STRIP_ALIGN_CELL):
    """Surfaces sol par faisceau sur maille 1 m (moyenne des points bas).

    Pour chaque faisceau et chaque maille : moyenne des points situés à moins
    de 0,5 m du minimum du faisceau dans la maille (robuste à la végétation
    résiduelle), mailles à ≥ 2 points seulement.

    Returns:
        (allcells, grid, key) : cellules triées, grid[faisceau, cellule] = Z
        moyen (NaN si absent) et clé de maille de CHAQUE point ; None si
        aucune surface n'est peuplée.
    """
    x0 = np.floor(np.min(x) / cell) * cell
    y0 = np.floor(np.min(y) / cell) * cell
    xi = ((x - x0) / cell).astype(np.int64)
    yi = ((y - y0) / cell).astype(np.int64)
    ny = int(yi.max()) + 1
    key = xi * ny + yi

    def _surface(k):
        m = inv == k
        kk, zz = key[m], z[m]
        order = np.argsort(kk, kind='stable')
        k_s, z_s = kk[order], zz[order]
        starts = np.flatnonzero(np.r_[True, k_s[1:] != k_s[:-1]])
        mins = np.minimum.reduceat(z_s, starts)
        ukey = k_s[starts]
        thr = mins[np.searchsorted(ukey, kk)]
        sel = np.flatnonzero(zz <= thr + 0.5)
        k2, z2 = kk[sel], zz[sel]
        cnt = np.bincount(k2)
        sums = np.bincount(k2, weights=z2)
        v = np.flatnonzero(cnt >= 2)
        return v, sums[v] / cnt[v]

    surfaces = [_surface(k) for k in range(n_sources)]
    populated = [c for c, _ in surfaces if len(c)]
    if not populated:
        return None
    allcells = np.unique(np.concatenate(populated))
    grid = np.full((n_sources, len(allcells)), np.nan)
    for k, (c, zs) in enumerate(surfaces):
        grid[k, np.searchsorted(allcells, c)] = zs
    return allcells, grid, key


def _strip_vertical_offsets(x, y, z, psid, cell=STRIP_ALIGN_CELL,
                            threshold=STRIP_ALIGN_THRESHOLD,
                            min_shared=STRIP_ALIGN_MIN_SHARED):
    """Mesure les offsets verticaux relatifs entre faisceaux d'une tuile.

    Méthode : surface sol par faisceau (moyenne des points à moins de 0,5 m du
    minimum de chaque maille de 1 m, robuste à la végétation résiduelle), puis
    offset de chaque faisceau = médiane de l'écart à la surface médiane de
    référence (itéré 3 fois), sur les seules cellules couvertes par au moins
    deux faisceaux. Ne corrige rien d'autre que le vertical : le calage
    horizontal mesuré sur les données est excellent (≤ 1 cm).

    Args:
        x, y, z, psid: coordonnées et PointSourceId des points sol.
        cell: taille de maille de comparaison (m).
        threshold: seuil (m) en dessous duquel un offset est ignoré.
        min_shared: cellules communes minimales pour valider un faisceau.

    Returns:
        dict {psid: offset} des offsets à SOUSTRAIRE (z - offset), ne
        contenant que les |offset| >= threshold ; vide si rien à corriger
        (faisceau unique, pas de recouvrement, tuile déjà alignée).
    """
    us, inv = np.unique(psid, return_inverse=True)
    if len(us) < 2:
        return {}
    built = _strip_surface_grid(x, y, z, inv, len(us), cell)
    if built is None:
        return {}
    allcells, grid, _key = built
    comparable = np.sum(~np.isnan(grid), axis=0) >= 2

    offsets = np.zeros(len(us))
    for _ in range(3):
        ref = np.nanmedian(grid + offsets[:, None], axis=0)
        for k in range(len(us)):
            m = ~np.isnan(grid[k]) & comparable
            if int(m.sum()) >= min_shared:
                offsets[k] = np.median(grid[k][m] - ref[m])

    return {int(p): round(float(offsets[k]), 3)
            for k, p in enumerate(us) if abs(offsets[k]) >= threshold}


def _rolling_median(values, window):
    """Médiane glissante (fenêtre tronquée aux bords) sur un petit vecteur."""
    n = len(values)
    if window <= 1 or n == 0:
        return np.array(values, dtype=np.float64)
    half = window // 2
    return np.array([np.median(values[max(0, i - half):i + half + 1])
                     for i in range(n)])


def _strip_jitter_offsets(x, y, z, psid, t, cell=STRIP_ALIGN_CELL,
                          bin_seconds=STRIP_JITTER_BIN,
                          smooth=STRIP_JITTER_SMOOTH,
                          min_cells=STRIP_JITTER_MIN_CELLS,
                          max_corr=STRIP_JITTER_MAX,
                          threshold=STRIP_ALIGN_THRESHOLD):
    """Mesure la gigue verticale intra-faisceau par fenêtres de temps GPS.

    Le calage constant retire un offset par faisceau, mais le décalage
    vertical peut aussi varier au fil d'une même passe (vibration capteur,
    bruit haute fréquence de la trajectoire). Chaque faisceau est découpé en
    fenêtres de temps GPS ; l'offset robuste de chaque fenêtre est mesuré
    contre la surface médiane des AUTRES faisceaux (maille 1 m, même
    sélection de points bas que le calage constant — en recouvrement à deux,
    une référence incluant le faisceau testé ne révélerait que la moitié du
    décalage), puis la série est lissée (médiane glissante) pour ne pas
    suivre le bruit de mesure. z doit déjà être corrigé des offsets constants.

    Args:
        x, y, z, psid: coordonnées, Z (constant-calé) et PointSourceId des
            points sol.
        t: temps GPS (s) de chaque point.
        cell: taille de maille de comparaison (m).
        bin_seconds: durée d'une fenêtre de temps (s).
        smooth: largeur de la médiane glissante (fenêtres).
        min_cells: cellules partagées minimales pour valider une fenêtre.
        max_corr: amplitude maximale d'une correction (garde-fou, m).
        threshold: amplitude de série sous laquelle un faisceau est
            considéré comme stable (pas de correction).

    Returns:
        dict {psid: (times, corrections)} des corrections à SOUSTRAIRE,
        interpolables linéairement au temps GPS de chaque point ; vide si
        rien de mesurable (faisceau unique, temps absent, pas de
        recouvrement).
    """
    us, inv = np.unique(psid, return_inverse=True)
    if len(us) < 2:
        return {}
    t = np.asarray(t, dtype=np.float64)
    if t.size != z.size or not np.isfinite(t).all():
        return {}
    built = _strip_surface_grid(x, y, z, inv, len(us), cell)
    if built is None:
        return {}
    allcells, grid, key = built
    covered = np.sum(~np.isnan(grid), axis=0) >= 2

    # Référence d'un faisceau = médiane des autres faisceaux.
    ref = np.full(grid.shape, np.nan)
    if len(us) == 2:
        ref[0], ref[1] = grid[1], grid[0]
    else:
        import warnings
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)  # tranches tout-NaN
            for k in range(len(us)):
                ref[k] = np.nanmedian(np.delete(grid, k, axis=0), axis=0)

    # Résidu vertical de chaque point contre la référence de son faisceau.
    idx = np.minimum(np.searchsorted(allcells, key), len(allcells) - 1)
    ref_point = ref[inv, idx]
    hit = (allcells[idx] == key) & covered[idx] & np.isfinite(ref_point)
    residual = np.full(np.asarray(z).shape, np.nan)
    residual[hit] = np.asarray(z)[hit] - ref_point[hit]

    # Origine de temps propre à chaque faisceau : des passes d'une même tuile
    # peuvent être espacées de plusieurs heures, les fenêtres restent ainsi
    # dense autour du vol réel (et un saut de semaine GPS ne crée pas de
    # géantes plages vides).
    origin = np.full(len(us), np.inf)
    np.minimum.at(origin, inv, t)
    tb = np.floor((t - origin[inv]) / bin_seconds).astype(np.int64)
    n_bins = int(tb.max()) + 1
    group = inv * n_bins + tb

    # Médiane robuste par groupe (faisceau, fenêtre) : les résidus finis sont
    # triés en tête de segment, les NaN (sans référence) sont ignorés.
    order = np.lexsort((np.where(np.isfinite(residual), residual, np.inf), group))
    g_s, r_s = group[order], residual[order]
    starts = np.flatnonzero(np.r_[True, g_s[1:] != g_s[:-1]])
    ends = np.r_[starts[1:], len(g_s)]

    gids, times_c, med = [], [], []
    for s, e in zip(starts, ends):
        finite = np.isfinite(r_s[s:e])
        if int(finite.sum()) < min_cells:
            continue
        gids.append(int(g_s[s]))
        med.append(float(np.median(r_s[s:e][finite])))
    if not gids:
        return {}

    gids = np.asarray(gids, dtype=np.int64)
    med = np.asarray(med, dtype=np.float64)
    result = {}
    for k in range(len(us)):
        selk = gids // n_bins == k
        if int(selk.sum()) < 3:  # série trop courte : correction non fiable
            continue
        b = gids[selk] % n_bins
        cs = np.clip(_rolling_median(med[selk], smooth), -max_corr, max_corr)
        if float(np.max(np.abs(cs))) < threshold:
            continue
        result[int(us[k])] = (origin[k] + (b + 0.5) * bin_seconds, cs)
    return result


def _apply_strip_jitter(psid, t, jitter):
    """Corrections de gigue interpolées au temps GPS de chaque point.

    Args:
        psid, t: PointSourceId et temps GPS (s) de chaque point.
        jitter: dict {psid: (times, corrections)} issu de _strip_jitter_offsets.

    Returns:
        array des corrections à SOUSTRAIRE (0 pour les points sans série).
    """
    corr = np.zeros(len(t))
    if not jitter:
        return corr
    psid = np.asarray(psid)
    t = np.asarray(t, dtype=np.float64)
    for p, (times, cs) in jitter.items():
        m = psid == p
        if m.any():
            corr[m] = np.interp(t[m], times, cs)
    return corr


def _strip_jitter_for_file(las_file, las, offsets):
    """Gigue temporelle d'un LAS sol (offsets constants mesurés), mémoïsée."""
    try:
        p = Path(las_file)
        cache_key = (str(p), p.stat().st_mtime_ns)
    except OSError:
        cache_key = (str(las_file), 0)
    if cache_key in _STRIP_JITTER_CACHE:
        return _STRIP_JITTER_CACHE[cache_key]
    jitter = {}
    try:
        t = np.asarray(las.gps_time, dtype=np.float64)
        psid = np.asarray(las.point_source_id)
    except AttributeError:
        t, psid = None, None  # dimensions absentes : pas de gigue mesurable
    if (t is not None and psid is not None
            and len(t) == len(las.points) == len(psid)):
        z = np.asarray(las.z, dtype=np.float64)
        if offsets:
            lut = np.zeros(65536)
            for p_, off in offsets.items():
                lut[int(p_) & 0xFFFF] = off
            z = z - lut[np.asarray(las.point_source_id, dtype=np.int64)]
        jitter = _strip_jitter_offsets(
            np.asarray(las.x, dtype=np.float64),
            np.asarray(las.y, dtype=np.float64),
            z, psid, t)
    _STRIP_JITTER_CACHE[cache_key] = jitter
    return jitter


def _strip_offsets_for_file(las_file, las):
    """Offsets de calage d'un LAS sol, mémoïsés par (chemin, mtime)."""
    try:
        p = Path(las_file)
        cache_key = (str(p), p.stat().st_mtime_ns)
    except OSError:
        cache_key = (str(las_file), 0)
    if cache_key in _STRIP_OFFSETS_CACHE:
        return _STRIP_OFFSETS_CACHE[cache_key]
    try:
        psid = np.asarray(las.point_source_id)
    except AttributeError:
        psid = None  # dimension absente (producteur tiers) : pas de calage
    if psid is not None and len(psid) == len(las.points):
        offsets = _strip_vertical_offsets(
            np.asarray(las.x, dtype=np.float64),
            np.asarray(las.y, dtype=np.float64),
            np.asarray(las.z, dtype=np.float64),
            psid)
    else:
        offsets = {}
    _STRIP_OFFSETS_CACHE[cache_key] = offsets
    return offsets


def _write_strip_align_sidecar(dtm_dir, basename, output_suffix, offsets,
                               jitter=None):
    """Consigne les calages appliqués (version, seuil, offsets et gigue).

    Le sidecar sert de suivi de cache : un DTM sans sidecar, ou produit avec
    une version/un seuil/des paramètres de gigue différents, est régénéré. Il
    est écrit même quand aucune correction n'a été appliquée, pour ne pas
    re-mesurer une tuile déjà connue comme bien alignée.
    """
    jitter_payload = {}
    for p, (times, cs) in (jitter or {}).items():
        cs = np.asarray(cs, dtype=np.float64)
        jitter_payload[str(p)] = {
            "bins": int(len(times)),
            "rms_m": round(float(np.sqrt(np.mean(cs ** 2))), 4),
            "max_m": round(float(np.max(np.abs(cs))), 4),
            "series_m": [round(float(c), 4) for c in cs],
        }
    payload = {
        "version": STRIP_ALIGN_VERSION,
        "threshold": STRIP_ALIGN_THRESHOLD,
        "offsets": offsets,
        "jitter_bin": STRIP_JITTER_BIN,
        "jitter_smooth": STRIP_JITTER_SMOOTH,
        "jitter": jitter_payload,
    }
    try:
        sidecar = Path(dtm_dir) / f"{basename}_dtm{output_suffix}_stripalign.json"
        sidecar.write_text(json.dumps(payload, ensure_ascii=False),
                           encoding="utf-8")
    except Exception as e:
        logger.warning(f"  Écriture sidecar calage faisceaux impossible: {e}")


def _strip_lidar_ext(path):
    """Extract base name from a LAZ/LAS file (mirrors pipeline._file_basename)."""
    name = Path(path).name
    for ext in ('.copc.laz', '.copc.las', '.laz', '.las'):
        if name.lower().endswith(ext):
            return name[:-len(ext)]
    return Path(path).stem


def parse_ign_classes(spec):
    """Convertit une liste de classes IGN (noms ou codes) en codes LAS triés.

    Args:
        spec: Chaîne séparée par virgules, ex. "sol,unclassified" ou "2,1".

    Returns:
        Liste triée de codes LAS uniques.

    Raises:
        ValueError: Si un élément n'est ni un nom connu ni un code LAS 0-255,
            ou si la liste est vide.
    """
    codes = set()
    for token in str(spec).split(","):
        token = token.strip().lower()
        if not token:
            continue
        if token in IGN_CLASS_NAMES:
            codes.add(IGN_CLASS_NAMES[token])
        else:
            try:
                code = int(token)
            except ValueError:
                raise ValueError(
                    f"Classe IGN inconnue: '{token}' "
                    f"(noms: {', '.join(sorted(IGN_CLASS_NAMES))} ou code LAS 0-255)")
            if not 0 <= code <= 255:
                raise ValueError(f"Code LAS hors bornes (0-255): {code}")
            codes.add(code)
    if not codes:
        raise ValueError("Aucune classe IGN fournie")
    return sorted(codes)


def ign_method_label(codes):
    """Étiquette de méthode encodant les classes IGN (ex. 'ign_1_2').

    'ign' seul = sol uniquement (code 2), rétrocompatible avec les fichiers de
    classification existants. Toute autre combinaison est encodée dans le nom
    pour invalider le cache et déclencher la reclassification.
    """
    codes = sorted(codes)
    if codes == [2]:
        return "ign"
    return "ign_" + "_".join(str(c) for c in codes)


def _create_ground_pipeline(input_laz, output_las, method, ign_codes=None):
    """Create a PDAL pipeline JSON for ground classification.

    All methods include a ReturnNumber/NumberOfReturns >= 1 filter to handle
    LiDAR HD files that may contain points with invalid return numbers.

    Pre-processing steps (PDAL recommended workflow):
    1. Reset Classification to 0
    2. ELM (Extended Local Minimum) — mark low outliers as noise (Classification=7)
    3. Statistical outlier removal
    4. Ground classification (SMRF or CSF)
    5. Extract ground points (Classification=2)

    Args:
        input_laz: Path to input LAZ/LAS file.
        output_las: Path to output classified LAS file.
        method: Ground classification method ('ign', 'smrf' or 'csf').
        ign_codes: LAS class codes to extract with the 'ign' method
            (default: [2] = sol). Multiple ranges on Classification are
            logically ORed by filters.range (documented PDAL semantics).

    Returns:
        JSON string of the PDAL pipeline.
    """
    # Common ReturnNumber filter for LiDAR HD compatibility
    return_filter = {
        "type": "filters.range",
        "limits": "ReturnNumber[1:],NumberOfReturns[1:]"
    }

    # Classification filter (ground points only)
    ground_filter = {
        "type": "filters.range",
        "limits": "Classification[2:2]"
    }

    # LiDAR HD IGN : le fichier est pré-classifié par le fournisseur.
    # On réutilise la classification telle quelle (mode pur) : les classes
    # extraites sont paramétrables — par défaut le sol seul (2), mais on peut
    # ajouter p.ex. unclassified (1) pour combler les trous sans retouche.
    # Les plages multiples sur Classification sont combinées en OU logique
    # par filters.range (sémantique PDAL documentée).
    if method == 'ign':
        codes = sorted(ign_codes) if ign_codes else [2]
        ign_filter = {
            "type": "filters.range",
            "limits": ",".join(f"Classification[{c}:{c}]" for c in codes)
        }
        pipeline = {
            "pipeline": [
                str(input_laz),
                return_filter,
                ign_filter,
                {
                    "type": "writers.las",
                    "filename": str(output_las),
                    "extra_dims": "all"
                }
            ]
        }
        return json.dumps(pipeline)

    # Reset Classification to 0 before preprocessing
    reset_classification = {
        "type": "filters.assign",
        "assignment": "Classification[:]=0"
    }

    # ELM (Extended Local Minimum) — mark low outliers as noise
    # Parameters tuned for rocky limestone terrain with low vegetation:
    # - cell=5.0m: fine resolution to capture rocky relief
    # - threshold=2.0m: high threshold to avoid marking rock outcrops as noise
    elm_filter = {
        "type": "filters.elm",
        "cell": 5.0,
        "threshold": 2.0
    }

    # Statistical outlier removal
    outlier_filter = {
        "type": "filters.outlier",
        "method": "statistical",
        "mean_k": 8,
        "multiplier": 3.0
    }

    # Method-specific ground classification filter
    if method == 'smrf':
        ground_step = {
            "type": "filters.smrf",
            "ignore": "Classification[7:7]",
            "slope": 1.0,
            "window": 16.0,
            "threshold": 0.5,
            "scalar": 1.25
        }
    elif method == 'csf':
        # resolution 1.0 m : un cloth à 0.5 m (= 4 M particules pour 1 km²)
        # rend la classification ~4× plus lente sans gain visible sur le MNT
        # (la résolution finale du MNT vient de la rasterisation, pas du cloth).
        ground_step = {
            "type": "filters.csf",
            "resolution": 1.0,
            "rigidness": 3,
            "smooth": True,
            "threshold": 0.5
        }
    else:
        raise ValueError(f"Méthode de classification inconnue: {method}")

    pipeline = {
        "pipeline": [
            str(input_laz),
            return_filter,
            reset_classification,
            elm_filter,
            outlier_filter,
            ground_step,
            ground_filter,
            {
                "type": "writers.las",
                "filename": str(output_las),
                "extra_dims": "all"
            }
        ]
    }
    return json.dumps(pipeline)


def create_smrf_pipeline(input_laz, output_las):
    """Create a PDAL pipeline JSON for SMRF ground classification."""
    return _create_ground_pipeline(input_laz, output_las, 'smrf')


def create_ign_pipeline(input_laz, output_las):
    """Create a PDAL pipeline JSON using the IGN supplier pre-classification."""
    return _create_ground_pipeline(input_laz, output_las, 'ign')


def create_csf_pipeline(input_laz, output_las):
    """Create a PDAL pipeline JSON for CSF ground classification."""
    return _create_ground_pipeline(input_laz, output_las, 'csf')


def validate_laz(laz_file):
    """Integrity check for a LAZ/LAS file.

    Verifies that both the header AND point data are readable. Some corrupted
    COPC files have valid headers but inaccessible point data (LazrsError:
    failed to fill whole buffer). Such files must be re-downloaded.

    Returns:
        True if file is readable and contains accessible points, False otherwise.
    """
    import laspy
    try:
        with laspy.open(str(laz_file)) as f:
            header = f.header
            point_count = header.point_count
            if point_count == 0:
                logger.error(f"  ✗ Fichier vide (0 points): {laz_file.name}")
                logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
                return False
            # Verify point data is actually accessible (not just header metadata)
            try:
                for _ in f.chunk_iterator(1000):
                    break  # Read just one chunk to confirm data integrity
            except Exception as chunk_err:
                logger.error(f"  ✗ Données inaccessibles (fichier corrompu?): {laz_file.name}")
                logger.error(f"    Erreur: {chunk_err}")
                logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
                return False
        return True
    except Exception:
        pass

    # Fallback: try PDAL (handles COPC v1.1 that laspy can't read)
    try:
        result = subprocess.run(
            ["pdal", "info", str(laz_file), "--summary"],
            capture_output=True, text=True, timeout=30
        )
        if result.returncode == 0:
            # Check point count from PDAL info output
            import json as _json
            try:
                info = _json.loads(result.stdout)
                count = info.get('summary', {}).get('num_points', 0)
                if count == 0:
                    logger.error(f"  ✗ Fichier vide (0 points PDAL): {laz_file.name}")
                    logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
                    return False
            except Exception:
                pass  # Can't parse — assume valid
            # Verify PDAL can actually read point data (not just header)
            try:
                test_result = subprocess.run(
                    ["pdal", "info", str(laz_file), "--point", "1"],
                    capture_output=True, text=True, timeout=60
                )
                if test_result.returncode != 0:
                    logger.error(f"  ✗ Données inaccessibles (PDAL): {laz_file.name}")
                    logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
                    return False
            except (subprocess.TimeoutExpired, FileNotFoundError):
                pass  # Timeout — assume valid, will fail later if corrupted
            return True
        logger.error(f"  ✗ Fichier illisible: {laz_file.name}")
        logger.error(f"    PDAL: {result.stderr.strip()[:200]}")
    except (subprocess.TimeoutExpired, FileNotFoundError):
        logger.error(f"  ✗ Impossible de vérifier le fichier: {laz_file.name}")
    logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
    return False


def _read_with_pdal(laz_file):
    """Read a LAZ/LAS file via PDAL when laspy fails (e.g. COPC v1.1).

    Returns a laspy.LasData object, or None on failure.
    """
    import subprocess
    import tempfile
    import os

    tmp_path = None
    try:
        # Convert COPC to LAS via PDAL, then read with laspy
        with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp:
            tmp_path = tmp.name

        pipeline = json.dumps({
            "pipeline": [
                str(laz_file),
                {"type": "writers.las", "filename": tmp_path}
            ]
        })

        result = subprocess.run(
            ["pdal", "pipeline", "--stdin"],
            input=pipeline, capture_output=True, text=True, timeout=300
        )

        if result.returncode != 0:
            logger.warning(f"  PDAL conversion échouée: {result.stderr[:200]}")
            return None

        import laspy
        las = laspy.read(tmp_path)
        if len(las.points) == 0:
            logger.warning(f"  PDAL: conversion réussie mais 0 points")
            return None
        return las

    except Exception as e:
        logger.warning(f"  PDAL fallback échoué: {e}")
        return None
    finally:
        # Toujours supprimer le LAS temporaire, même en cas d'exception
        # (timeout subprocess, laspy.read échoué…) : sinon fuite de ~Go.
        if tmp_path:
            try:
                os.unlink(tmp_path)
            except Exception:
                pass


def detect_ground_method(laz_file):
    """Detect the best ground classification method based on point cloud statistics.

    Auto-selects between SMRF and CSF:
    - SMRF: fast, robust for most natural terrain (PDAL recommended default)
    - CSF: cloth simulation, better for complex/urban terrain

    Falls back to SMRF if the file cannot be read or attributes are missing.

    Args:
        laz_file: Path to input LAZ/LAS file.

    Returns:
        String: 'ign', 'smrf' or 'csf'
    """
    import laspy

    # Try laspy first, then PDAL for COPC files
    las = None
    try:
        las = laspy.read(str(laz_file))
    except Exception as e:
        logger.warning(f"  laspy: {e}")
        logger.info(f"  → Lecture via PDAL pour auto-détection...")
        las = _read_with_pdal(laz_file)

    if las is None:
        logger.info(f"  → Méthode: SMRF (défaut — lecture impossible)")
        return 'smrf'

    total_points = len(las.points)
    if total_points == 0:
        logger.warning(f"  Nuage vide (0 points) — méthode par défaut: SMRF")
        return 'smrf'

    # LiDAR HD IGN : les données livrées sont pré-classifiées par le fournisseur
    # (classe 2 = sol). C'est la base la plus rapide (~10 s) et de référence.
    # Le MNT est ensuite complété par le retour le plus bas par cellule +
    # interpolation (voir create_dtm_fast), ce qui « rattrape » les trous de la
    # pré-classification (forêt dense / relief). On la préfère donc dès qu'une
    # part raisonnable des points est classée sol, plutôt que de refiltrer.
    try:
        cls = np.asarray(las.classification, dtype=np.int32)
        ground_ratio = float(np.mean(cls == 2))
    except Exception:
        ground_ratio = 0.0
    if ground_ratio >= 0.2:
        logger.info(f"  → Méthode: IGN (pré-classification fournisseur — "
                    f"{ground_ratio * 100:.1f}% de points classe 2)")
        return 'ign'

    z = np.array(las.z)

    # Height variance (always available)
    z_std = float(np.std(z))
    z_range = float(np.max(z) - np.min(z))

    # Try to get NumberOfReturns (may not exist in all point formats)
    single_return_ratio = 0.0
    try:
        num_returns = np.array(las.NumberOfReturns)
        single_return_count = int(np.sum(num_returns == 1))
        single_return_ratio = single_return_count / total_points if total_points > 0 else 0
    except AttributeError:
        logger.debug("  NumberOfReturns non disponible — utilisation de la variance Z uniquement")

    logger.info(f"  Analyse du nuage: {total_points:,} points, "
                f"ratio_retours_uniques={single_return_ratio:.2f}, "
                f"écart_Z={z_std:.1f}m, amplitude_Z={z_range:.1f}m")

    # Decision logic:
    # - High single-return ratio (>0.6) → urban (buildings, roads) → CSF (cloth simulation)
    # - High elevation variance (>30m) → complex/mountainous terrain → CSF
    # - Default → SMRF (fast, robust for most natural terrain)
    if single_return_ratio > 0.6:
        method = 'csf'
        reason = f"ratio retours uniques={single_return_ratio:.2f} > 0.6 → milieu urbain"
    elif z_std > 30:
        method = 'csf'
        reason = f"écart_Z={z_std:.1f}m > 30m → terrain complexe"
    else:
        method = 'smrf'
        reason = f"terrain naturel standard"

    logger.info(f"  → Méthode: {method.upper()} ({reason})")
    return method


def classify_ground(laz_file, temp_dir, method='auto', force=False, ign_classes="sol"):
    """Classify ground points using PDAL ground classification filter.

    Args:
        laz_file: Path to input LAZ/LAS file.
        temp_dir: Directory for temporary files (pipeline.json, ground.las).
        method: Ground classification method ('auto', 'ign', 'smrf' or 'csf').
        force: If True, reclassify even if output file already exists.
        ign_classes: Classes LAS extraites par la méthode IGN (noms ou codes
            séparés par virgules, ex. "sol,unclassified"). Ignoré pour les
            autres méthodes.

    Returns:
        Path to classified ground LAS file, or None on failure.
    """
    import laspy

    # Auto-detect method if requested
    if method == 'auto':
        method = detect_ground_method(laz_file)
        logger.info(f"  Classification sol: {method.upper()} (auto)")
    else:
        logger.info(f"  Classification sol: {method.upper()} (forcé)")

    # Les classes IGN sont encodées dans le nom de fichier (ex. ign_1_2)
    # pour qu'un changement de classes invalide le cache et déclenche la
    # reclassification.
    ign_codes = parse_ign_classes(ign_classes) if method == 'ign' else None
    method_label = ign_method_label(ign_codes) if ign_codes else method

    laz_base = _strip_lidar_ext(laz_file)

    output_las = temp_dir / f"{laz_base}_ground_{method_label}.las"

    if output_las.exists() and not force:
        logger.info(f"  Classification {method.upper()} déjà effectuée — fichier existant réutilisé")
        return output_las

    if force and output_las.exists():
        logger.info(f"  Reclassification forcée — suppression de {output_las.name}")
        output_las.unlink()

    pipeline_json = _create_ground_pipeline(laz_file, output_las, method, ign_codes=ign_codes)
    pipeline_file = temp_dir / f"pipeline_{method_label}.json"

    with open(pipeline_file, 'w') as f:
        f.write(pipeline_json)

    try:
        subprocess.run(
            ["pdal", "pipeline", str(pipeline_file)],
            capture_output=True, check=True
        )
        # Verify that ground file has points
        if output_las.exists() and output_las.stat().st_size < 100:
            logger.error(f"  ✗ Fichier ground vide (taille < 100 octets)")
            output_las.unlink(missing_ok=True)
            # Fallback: si la méthode ne produit aucun point sol, réessayer avec SMRF
            if method in ('csf', 'ign'):
                return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)
            return None
        logger.info(f"  ✓ Classification sol {method.upper()} terminée")
        return output_las
    except subprocess.CalledProcessError as e:
        error_msg = e.stderr.decode() if e.stderr else str(e)
        logger.warning(f"  ✗ Erreur classification PDAL ({method.upper()}): {error_msg}")

        # Fallback: si CSF ou la pré-classification échouent, réessayer avec SMRF
        if method in ('csf', 'ign'):
            return _fallback_to_smrf(laz_file, temp_dir, laz_base, force, source=method_label)

        # Try repairing file with laspy if PDAL fails on EVLR/VLR
        if 'VLR' in error_msg or 'Invalid' in error_msg:
            logger.info(f"  → Tentative de réparation du fichier avec laspy...")
            repaired_las = temp_dir / f"{laz_base}_repaired.las"
            if _repair_laz_with_laspy(laz_file, repaired_las):
                # Retry PDAL pipeline with repaired file
                pipeline_json = _create_ground_pipeline(repaired_las, output_las, method)
                with open(pipeline_file, 'w') as f:
                    f.write(pipeline_json)
                try:
                    subprocess.run(
                        ["pdal", "pipeline", str(pipeline_file)],
                        capture_output=True, check=True
                    )
                    logger.info(f"  ✓ Classification sol {method.upper()} terminée (fichier réparé)")
                    return output_las
                except subprocess.CalledProcessError as e2:
                    error_msg2 = e2.stderr.decode() if e2.stderr else str(e2)
                    logger.error(f"  ✗ Échec classification même après réparation: {error_msg2}")
            else:
                logger.error(f"  ✗ Impossible de réparer le fichier")
        return None


def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False, source='csf'):
    """Retry ground classification with SMRF when CSF/IGN fails.

    CSF (Cloth Simulation Filter) can fail on certain terrain types where
    SMRF (Simple Morphological Filter) succeeds, and a file without usable
    pre-classification produces an empty ground extract. This fallback ensures
    processing continues even when the selected method fails.

    Args:
        laz_file: Path to input LAZ/LAS file.
        temp_dir: Directory for temporary files.
        laz_base: Base name for the file.
        force: If True, reclassify even if output exists.
        source: Method that failed ('csf' or 'ign').

    Returns:
        Path to classified ground LAS file, or None on failure.
    """
    logger.info(f"  → Basculement {source.upper()} → SMRF (fallback)")

    # Clean up failed output if it exists
    failed_output = temp_dir / f"{laz_base}_ground_{source}.las"
    if failed_output.exists():
        failed_output.unlink(missing_ok=True)

    output_las = temp_dir / f"{laz_base}_ground_smrf.las"

    if output_las.exists() and not force:
        logger.info(f"  Classification SMRF déjà existante — fichier réutilisé")
        return output_las

    pipeline_json = _create_ground_pipeline(laz_file, output_las, 'smrf')
    pipeline_file = temp_dir / "pipeline_smrf.json"

    with open(pipeline_file, 'w') as f:
        f.write(pipeline_json)

    try:
        subprocess.run(
            ["pdal", "pipeline", str(pipeline_file)],
            capture_output=True, check=True
        )
        if output_las.exists() and output_las.stat().st_size < 100:
            logger.error(f"  ✗ Fichier ground SMRF vide (taille < 100 octets)")
            output_las.unlink(missing_ok=True)
            return None
        logger.info(f"  ✓ Classification sol SMRF terminée (fallback depuis CSF)")
        return output_las
    except subprocess.CalledProcessError as e:
        error_msg = e.stderr.decode() if e.stderr else str(e)
        logger.error(f"  ✗ Échec classification SMRF (fallback): {error_msg}")
        return None


def _repair_laz_with_laspy(input_laz, output_las):
    """Try to repair a corrupt LAZ file by re-reading with laspy and saving as LAS.

    Works around PDAL errors like 'Invalid Extended VLR size' by stripping
    problematic VLR/EVLR metadata during re-save.

    Args:
        input_laz: Path to corrupt LAZ/LAS file.
        output_las: Path for repaired LAS output.

    Returns:
        True if repair succeeded, False otherwise.
    """
    import laspy
    try:
        las = laspy.read(str(input_laz))
        las.write(str(output_las))
        logger.info(f"  ✓ Fichier réparé via laspy ({len(las.points):,} points)")
        return True
    except Exception as e:
        logger.warning(f"  ✗ Réparation laspy échouée: {e}")
        return False


def _interpolate_holes(dtm, downsample=8):
    """Fill remaining NaN holes with a terrain-aware surface interpolation.

    In complex / rocky terrain the ground under-classification leaves interior
    holes far too large for a 1 m gap fill, which otherwise become flat
    nearest-neighbor patches in the downstream layers. This helper triangulates
    the valid cells on a downsampled grid (linear, nearest as a fallback for
    cells outside the data hull) and bilinearly upsamples the result, keeping
    the operation fast even for large rasters.

    Args:
        dtm: 2-D float array (may contain NaN holes).
        downsample: Coarsening factor for the interpolation grid.

    Returns:
        Tuple (filled_array, filled_count).
    """
    holes = np.isnan(dtm)
    if not holes.any():
        return dtm, 0
    valid = ~holes
    if not valid.any():
        return dtm, 0

    height, width = dtm.shape
    step = max(1, downsample)
    coarse = dtm[::step, ::step].astype(np.float64)
    c_valid = ~np.isnan(coarse)
    c_holes = np.isnan(coarse)
    if not c_holes.any() or not c_valid.any():
        return dtm, 0

    from scipy.interpolate import griddata
    from scipy.ndimage import map_coordinates

    cy, cx = np.where(c_valid)
    c_coords = np.column_stack([cx, cy]).astype(np.float64)
    c_vals = coarse[c_valid]
    hy, hx = np.where(c_holes)
    h_coords = np.column_stack([hx, hy]).astype(np.float64)

    interp = griddata(c_coords, c_vals, h_coords, method='linear')
    bad = np.isnan(interp)
    if bad.any():
        interp[bad] = griddata(c_coords, c_vals, h_coords[bad], method='nearest')
    coarse_filled = coarse.copy()
    coarse_filled[c_holes] = interp

    # Coarse cell i represents fine column/row i*step, so fine index c maps to
    # coarse coordinate c/step (no half-cell offset).
    rows = np.arange(height) / step
    cols = np.arange(width) / step
    grid_y, grid_x = np.meshgrid(rows, cols, indexing='ij')
    upsampled = map_coordinates(coarse_filled, [grid_y, grid_x], order=1)

    filled = dtm.copy()
    filled[holes] = upsampled[holes]
    return filled, int(holes.sum())


# ------------------------------------------------------------
# Raccord des bords (edge buffer avec les tuiles adjacentes)
# ------------------------------------------------------------
# Les visualisations à grand noyau (openness/SVF : rayons jusqu'à 100 m,
# LRM : 15 m) tronquent leur fenêtre au bord de la dalle : les rendus
# présentent alors une bande d'artefacts à chaque changement de tuile. Avec
# edge_buffer > 0, le MNT est rastérisé sur la dalle nominale 1 km ÉTENDUE
# d'une bande de `edge_buffer` mètres remplie avec les points sol des 8
# tuiles LAZ voisines (lecture PDAL en flux : découpe + filtre de classes).
# Les visualisations voient ainsi le terrain réel au-delà du bord, et les
# images finales sont recadrées sur la dalle exacte (cf. rendering.py).
EDGE_BUFFER_TAG = "LIDAR_EDGE_BUFFER"  # tag GeoTIFF : tampon utilisé (m)

# Sous-dossier de input/ où sont isolées les dalles voisines téléchargées
# pour le seul raccord des bords : elles ne font pas partie du corpus de
# tuiles à rendre (les scans globaux n'énumèrent que input/ à plat).
EDGE_NEIGHBORS_DIRNAME = "edge_neighbors"

# Décalages des 8 voisines d'une dalle (col, row) en km
_NEIGHBOR_OFFSETS = [(-1, -1), (0, -1), (1, -1), (-1, 0),

                     (1, 0), (-1, 1), (0, 1), (1, 1)]


def _tile_coords(name):
    """Coordonnées (col, row) en km d'un nom de fichier LHD, ou None."""
    from .index import parse_basename_coords
    return parse_basename_coords(Path(name).name)


def _neighbor_laz_files(source_laz):
    """Liste les 8 fichiers LAZ/LAS adjacents à `source_laz`.

    Recherche dans le dossier de la source, puis dans le sous-dossier
    edge_neighbors/ : les dalles voisines téléchargées pour le seul raccord
    des bords y sont isolées — sinon elles s'accumulent à plat dans input/ et
    chaque passe globale annexe un anneau de tuiles à rendre en plus.
    """
    coords = _tile_coords(source_laz)
    if coords is None:
        return []
    col, row = coords
    directory = Path(source_laz).parent
    search_dirs = [directory, directory / EDGE_NEIGHBORS_DIRNAME]
    neighbors = []
    for dcol, drow in _NEIGHBOR_OFFSETS:
        nc, nr = col + dcol, row + drow
        # Noms LHD : col/row en km sur 4 chiffres complétés (ex. 0638_6628) ;
        # on accepte aussi la variante sans remplissage pour les dalles exotiques.
        names = {f"LHD_FXX_{nc:04d}_{nr:04d}", f"LHD_FXX_{nc}_{nr}"}
        matches = []
        for search_dir in search_dirs:
            for name in names:
                matches += list(search_dir.glob(f"{name}_*.las"))
                matches += list(search_dir.glob(f"{name}_*.laz"))
        if matches:
            neighbors.append(sorted(matches)[0])
        else:
            logger.debug(f"    Voisine absente : LHD_FXX_{nc:04d}_{nr:04d} (bande de bord vide)")
    return neighbors


def _neighbor_ground_points(source_laz, bounds, classes):
    """Points sol des tuiles voisines dans `bounds` (raccord des bords).

    Lecture PDAL en flux par voisine : découpe sur l'emprise étendue puis
    filtre de classes (mêmes codes que le MNT). Best-effort : une voisine
    illisible ou absente est ignorée — la bande correspondante reste vide.

    Args:
        source_laz: LAZ de la tuile traitée (repère pour trouver les voisines).
        bounds: (min_x, min_y, max_x, max_y) de l'emprise étendue.
        classes: codes LAS à extraire (ex. [2] = sol).

    Returns:
        (xs, ys, zs) concaténés (tableaux vides si aucune voisine).
    """
    import tempfile

    neighbors = _neighbor_laz_files(source_laz)
    if not neighbors:
        return (np.empty(0),) * 3

    min_x, min_y, max_x, max_y = bounds
    codes = sorted(set(int(c) for c in classes)) or [2]
    limits = ",".join(f"Classification[{c}:{c}]" for c in codes)

    xs, ys, zs = [], [], []
    found = 0
    for neighbor in neighbors:
        tmp_path = None
        try:
            with tempfile.NamedTemporaryFile(suffix='.las', delete=False) as tmp:
                tmp_path = tmp.name
            pipeline = json.dumps({
                "pipeline": [
                    {"type": "readers.las", "filename": str(neighbor)},
                    {"type": "filters.crop",
                     "bounds": f"([{min_x},{max_x}],[{min_y},{max_y}])"},
                    {"type": "filters.range", "limits": limits},
                    {"type": "writers.las", "filename": tmp_path},
                ]
            })
            result = subprocess.run(
                ["pdal", "pipeline", "--stdin"],
                input=pipeline, capture_output=True, text=True, timeout=300
            )
            if result.returncode != 0:
                raise RuntimeError(result.stderr[:200])
            import laspy
            las = laspy.read(tmp_path)
            if len(las.points) > 0:
                xs.append(np.asarray(las.x, dtype=np.float64))
                ys.append(np.asarray(las.y, dtype=np.float64))
                zs.append(np.asarray(las.z, dtype=np.float64))
                found += 1
        except Exception as e:
            logger.debug(f"    Voisine {Path(neighbor).name} ignorée : {e}")
        finally:
            if tmp_path:
                try:
                    Path(tmp_path).unlink(missing_ok=True)
                except Exception:
                    pass

    if not found:
        logger.warning("    Aucune voisine lisible — bande de bord vide")
        return (np.empty(0),) * 3
    logger.info(f"    Raccord bords : {found} voisine(s), "
                f"{sum(len(a) for a in xs):,} pts sol")
    return np.concatenate(xs), np.concatenate(ys), np.concatenate(zs)


def read_dtm_edge_buffer(dtm_path):
    """Tampon de raccord enregistré dans un DTM (m ; 0 si absent/illisible)."""
    try:
        with rasterio.open(dtm_path) as src:
            tags = src.tags()
        return float(tags.get(EDGE_BUFFER_TAG, 0.0) or 0.0)
    except Exception:
        return 0.0


def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
                    output_suffix="", source_laz=None,
                    pure=False, strip_align=True, edge_buffer=0.0,
                    neighbor_classes=None):
    """Create DTM using fast binning method with gap filling.

    Args:
        las_file: Path to classified ground LAS file.
        basename: Base name for output file.
        dtm_dir: Directory for output DTM GeoTIFF.
        resolution: Grid resolution in meters per pixel.
        force: If True, regenerate even if DTM already exists.
        output_suffix: Suffix for output filename (e.g. '_r0p2' for additional resolutions).
        source_laz: Optionnel : LAZ complet de la tuile traitée (repère pour
            lire les 8 tuiles voisines du raccord de bords, cf. edge_buffer).
        pure: Sans effet (conservé pour compatibilité). Fonctionnement
            historique rétabli : petits trous comblés par fillnodata, grands
            trous laissés en nodata (rendus en noir dans les rendus).
        strip_align: Si True (défaut), mesure et corrige les écarts verticaux
            entre faisceaux de vol (PointSourceId) avant rastérisation : offsets
            constants ≥ STRIP_ALIGN_THRESHOLD (0,5 cm) puis gigue intra-faisceau
            par fenêtres de temps GPS ; tout est consigné dans un sidecar
            *_dtm*_stripalign.json.
        edge_buffer: Raccord des bords en mètres (0 = désactivé). Le MNT couvre
            alors la dalle nominale 1 km étendue de cette bande, remplie avec
            les points sol des 8 tuiles LAZ voisines (source_laz requis) ; les
            visualisations calculent sur l'emprise étendue puis les images sont
            recadrées sur la dalle exacte (rendering.py). Le tampon est inscrit
            dans le tag GeoTIFF LIDAR_EDGE_BUFFER pour l'invalidation du cache.
        neighbor_classes: codes LAS extraits chez les voisines (défaut : les
            classes IGN du MNT, ex. [2]). Les voisines sont lues dans leur
            pré-classification fournisseur, même si la tuile centrale est
            classée SMRF/CSF (bande de contexte, quelques cm d'écart au pire).

    Returns:
        Path to output DTM GeoTIFF, or None on failure.
    """
    output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"

    if output_tif.exists() and not force:
        logger.info(f"  DTM déjà existant — fichier réutilisé: {output_tif.name}")
        return output_tif

    import laspy

    logger.info("  → Génération DTM...")

    try:
        las = laspy.read(str(las_file))
    except Exception as e:
        # laspy can't read COPC v1.1 — try PDAL conversion
        logger.warning(f"  laspy: {e}")
        logger.info(f"  → Conversion via PDAL pour lecture COPC...")
        las = _read_with_pdal(las_file)
        if las is None:
            logger.error(f"  ✗ Impossible de lire {las_file.name}")
            return None

    if len(las.points) == 0:
        logger.error(f"  ✗ Fichier vide (0 points): {las_file.name}")
        return None

    # Calage vertical des faisceaux de vol avant rastérisation : best-effort,
    # en cas d'échec de la mesure on continue non calé (jamais d'abort).
    strip_offsets = {}
    strip_jitter = {}
    gps_time = None
    if strip_align:
        try:
            strip_offsets = _strip_offsets_for_file(las_file, las)
        except Exception as e:
            logger.warning(f"  Mesure du calage faisceaux impossible ({e}) — MNT non calé")
            strip_offsets = {}
        if strip_offsets:
            logger.info("    Calage faisceaux : " + ", ".join(
                f"PSID {p} {off:+.3f} m" for p, off in sorted(strip_offsets.items())))
        else:
            logger.debug("    Calage faisceaux : aucun écart >= "
                         f"{STRIP_ALIGN_THRESHOLD * 100:.1f} cm, rien à corriger")
        # 2ᵉ passe : gigue intra-faisceau (temps GPS requis, ignorée sinon).
        try:
            gps_time = np.asarray(las.gps_time, dtype=np.float64)
        except AttributeError:
            gps_time = None
        if gps_time is not None and len(gps_time) == len(las.points):
            try:
                strip_jitter = _strip_jitter_for_file(las_file, las, strip_offsets)
            except Exception as e:
                logger.warning(f"  Mesure de la gigue intra-faisceau impossible ({e}) — gigue non corrigée")
                strip_jitter = {}
        if strip_jitter:
            for p_, (jt, jc) in sorted(strip_jitter.items()):
                logger.info(f"    Gigue PSID {p_} : ±{np.max(np.abs(jc)) * 100:.1f} cm "
                            f"(rms {np.sqrt(np.mean(np.asarray(jc) ** 2)) * 100:.1f} cm, "
                            f"{len(jt)} fenêtres de {STRIP_JITTER_BIN:g} s)")
        else:
            logger.debug("    Gigue intra-faisceau : rien à corriger")

    try:

        min_x, max_x = float(las.header.min[0]), float(las.header.max[0])
        min_y, max_y = float(las.header.min[1]), float(las.header.max[1])

        width = int(np.ceil((max_x - min_x) / resolution))
        height = int(np.ceil((max_y - min_y) / resolution))

        # Raccord des bords : emprise = dalle nominale 1 km (alignée sur la
        # grille multi-tuiles) + bande de edge_buffer mètres remplie par les
        # points sol des voisines. Sinon : bornes de l'en-tête (historique).
        used_edge_buffer = 0.0
        if edge_buffer > 0:
            coords = _tile_coords(source_laz or las_file)
            if coords is not None:
                buffer_px = max(1, int(round(edge_buffer / resolution)))
                buffer_m = buffer_px * resolution
                col_km, row_km = coords
                # Grille LHD : (col, row) = coin nord-ouest en km →
                # X ∈ [col, col+1] km, Y ∈ [row-1, row] km (bord nord = row).
                min_x = float(col_km) * 1000.0
                max_x = min_x + 1000.0
                max_y = float(row_km) * 1000.0
                min_y = max_y - 1000.0
                ext_bounds = (min_x - buffer_m, min_y - buffer_m,
                              max_x + buffer_m, max_y + buffer_m)
                width = int(round(1000.0 / resolution)) + 2 * buffer_px
                height = width
                min_x, min_y, max_x, max_y = ext_bounds
                used_edge_buffer = float(edge_buffer)
            else:
                logger.warning("    Raccord des bords impossible : nom de "
                               f"fichier non LHD ({basename}) — tuile seule")

        logger.debug(f"    Bounds: X[{min_x:.1f}, {max_x:.1f}] Y[{min_y:.1f}, {max_y:.1f}]")
        logger.debug(f"    Grid: {width}x{height} pixels ({len(las.points):,} points)")
        logger.info(f"    Rasterisation {width}x{height} ({len(las.points):,} points)...")

        xs = np.asarray(las.x, dtype=np.float64)
        ys = np.asarray(las.y, dtype=np.float64)
        zs = np.asarray(las.z, dtype=np.float64)
        if strip_offsets:
            lut = np.zeros(65536)
            for p, off in strip_offsets.items():
                lut[int(p) & 0xFFFF] = off
            zs = zs - lut[np.asarray(las.point_source_id, dtype=np.int64)]
        if strip_jitter and gps_time is not None:
            zs = zs - _apply_strip_jitter(las.point_source_id, gps_time, strip_jitter)

        # Points sol des tuiles voisines dans la bande de raccord (best-effort,
        # non calés par faisceau : bande de contexte, l'image finale est
        # recadrée sur la dalle avant livraison).
        if used_edge_buffer > 0 and source_laz is not None:
            nx, ny, nz = _neighbor_ground_points(
                source_laz, (min_x, min_y, max_x, max_y),
                neighbor_classes if neighbor_classes is not None else [2])
            if len(nx):
                xs = np.concatenate([xs, nx])
                ys = np.concatenate([ys, ny])
                zs = np.concatenate([zs, nz])

        stat = binned_statistic_2d(
            xs, ys, zs,
            statistic='mean',
            bins=[width, height],
            range=[[min_x, max_x], [min_y, max_y]]
        )

        dtm = stat.statistic.T
        dtm = dtm[::-1, :]  # Flip Y so north is at top

        # Comblement « historique » (fonctionnement d'origine, rétabli) :
        # seuls les petits trous proches des données sont remplis ; les grands
        # trous restent en nodata et apparaissent en noir dans les rendus.
        # Volontairement PAS de plancher au retour le plus bas : sous canopée
        # dense ce retour est la végétation, qui imprimerait les arbres dans
        # le MNT.

        # Fill small gaps (< 1 m from data) precisely — comme avant
        nan_count = np.count_nonzero(np.isnan(dtm))
        if nan_count > 0:
            total = dtm.size
            nan_pct = 100.0 * nan_count / total
            logger.info(f"    {nan_count:,} pixels sans données ({nan_pct:.1f}%)")

            max_gap_pixels = max(1, int(1.0 / resolution))
            from rasterio.fill import fillnodata
            valid_mask = ~np.isnan(dtm)
            dtm_filled = fillnodata(dtm, mask=valid_mask, max_search_distance=max_gap_pixels)
            small_gap_mask = np.isnan(dtm) & ~np.isnan(dtm_filled)
            filled_count = np.count_nonzero(small_gap_mask)
            if filled_count > 0:
                dtm = np.where(small_gap_mask, dtm_filled, dtm)
                logger.info(f"    {filled_count:,} petits trous comblés (< {max_gap_pixels}px)")

        # Save as GeoTIFF
        output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"
        transform = from_bounds(min_x, min_y, max_x, max_y, width, height)

        with rasterio.open(
            output_tif, 'w',
            driver='GTiff', height=height, width=width,
            count=1, dtype='float32',
            crs='EPSG:2154', transform=transform,
            nodata=float('nan'),
            compress='lzw'
        ) as dst:
            dst.write(dtm.astype('float32'), 1)
            if used_edge_buffer > 0:
                # Tampon de raccord inscrit dans le fichier : changement de
                # --edge-buffer ⇒ invalidation automatique du cache DTM.
                dst.update_tags(**{EDGE_BUFFER_TAG: f"{used_edge_buffer:g}"})

        if strip_align:
            _write_strip_align_sidecar(dtm_dir, basename, output_suffix,
                                       strip_offsets, strip_jitter)
        logger.info(f"    ✓ DTM créé: {output_tif.name}")
        return output_tif

    except Exception as e:
        logger.error(f"  ✗ Erreur DTM: {e}", exc_info=True)
        return None