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


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  # noqa: ensure available

    # 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())


def _min_return_grid(laz_file, width, height, bounds, chunk_size=2_000_000):
    """Rasterize the per-cell minimum z (lowest return) of the full point cloud.

    In complex/forested terrain the ground is under-classified, leaving DTM
    holes. Filling them with the *lowest measured return* of the cell (Wack &
    Wimmer 2002) recovers a real ground surface (forest floor, rock, clearing)
    instead of a pure interpolation. The read is streamed in chunks so memory
    stays bounded to the output grid regardless of the point count.

    Args:
        laz_file: Path to the full (unclassified) LAZ/LAS file.
        width, height: Output grid dimensions (pixels).
        bounds: (min_x, min_y, max_x, max_y) the grid covers.
        chunk_size: Points per streaming chunk.

    Returns:
        (height, width) float32 array of per-cell min z (NaN where no point).
    """
    import laspy
    min_x, min_y, max_x, max_y = bounds
    grid = np.full((height, width), np.nan, dtype=np.float32)
    rng = [[min_x, max_x], [min_y, max_y]]

    def process(points):
        if len(points) == 0:
            return
        x = np.asarray(points.x, dtype=np.float64)
        y = np.asarray(points.y, dtype=np.float64)
        z = np.asarray(points.z, dtype=np.float64)
        st = binned_statistic_2d(x, y, z, statistic='min',
                                 bins=[width, height], range=rng)
        # Match the DTM convention: .T then flip Y (north at top).
        cell_min = st.statistic.T[::-1, :].astype(np.float32)
        # fmin ignores NaN so cells without a point in this chunk stay NaN.
        np.fmin(grid, cell_min, out=grid)

    try:
        with laspy.open(str(laz_file)) as las:
            for chunk in las.chunk_iterator(chunk_size):
                process(chunk)
    except Exception as e:
        logger.warning(f"    Lecture streaming impossible ({e}) — lecture complète")
        las = _read_with_pdal(laz_file)
        if las is None:
            return grid
        process(las)
    return grid


def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False,
                    output_suffix="", source_laz=None, bare_earth=False,
                    pure=False):
    """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 : chemin du LAZ complet (non classé). Utilisé
            uniquement avec bare_earth (plancher au retour le plus bas).
        bare_earth: If True, pull the DTM down to the lowest measured return of
            each cell (bare-earth floor). This requalifies the lowest point of
            every column as terrain, recovering the ground under dense
            vegetation / steep relief that the ground classifier rejected.
        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).

    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

    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))

        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)...")

        stat = binned_statistic_2d(
            las.x, las.y, las.z,
            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.
        # Le plancher au retour le plus bas n'est appliqué qu'à la demande
        # explicite (--bare-earth).
        if bare_earth and source_laz is not None:
            min_grid = _min_return_grid(source_laz, width, height,
                                        (min_x, min_y, max_x, max_y))
            # Cellules sans sol mesuré (NaN) : le retour le plus bas devient
            # la mesure — sinon la comparaison NaN est fausse et la cellule
            # retombe sur l'interpolation en fin de passe.
            has_min = ~np.isnan(min_grid)
            lower = has_min & (min_grid < dtm)
            lower |= has_min & np.isnan(dtm)
            dtm = np.where(lower, min_grid, dtm)
            logger.info(f"    Sol nu : {int(lower.sum()):,} cellules raménées au retour le plus bas")

        # 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)

        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