"""Terrain visualization functions for LiDAR archaeological analysis.

Each function takes (dem_file, basename, vis_dir, resolution) as explicit
parameters and returns the path to the output GeoTIFF file, or None on error.
"""

import logging
import time
from pathlib import Path

import numpy as np
import rasterio
from scipy import ndimage
from scipy.ndimage import generic_filter, gaussian_filter, uniform_filter, minimum_filter
from scipy.stats import binned_statistic_2d

from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter

logger = logging.getLogger("lidar")

# Use CuPy array module when available
if HAS_GPU:
    import cupy as cp
    xp = cp
else:
    xp = np


def _save_tif(output_path, data, transform, crs, dtype='float32', count=1):
    """Helper to save a 2D or 3D array as GeoTIFF."""
    if data.ndim == 2:
        height, width = data.shape
        with rasterio.open(
            output_path, 'w', driver='GTiff',
            height=height, width=width, count=count,
            dtype=dtype, crs=crs, transform=transform, compress='lzw'
        ) as dst:
            dst.write(data.astype(dtype), 1)
    elif data.ndim == 3:
        bands, height, width = data.shape
        with rasterio.open(
            output_path, 'w', driver='GTiff',
            height=height, width=width, count=bands,
            dtype=dtype, crs=crs, transform=transform, compress='lzw'
        ) as dst:
            for i in range(bands):
                dst.write(data[i].astype(dtype), i + 1)


def _read_dem(dem_file):
    """Read DEM file and return (data, transform, crs)."""
    with rasterio.open(dem_file) as src:
        return src.read(1), src.transform, src.crs


# ============================================================
# Core terrain visualizations
# ============================================================

def generate_hillshade(dem_file, basename, vis_dir, resolution):
    """Generate multi-directional hillshade (NW, NE, SW, SE)."""
    logger.info("  → Hillshade multidirectionnel...")
    t0 = time.time()
    output = vis_dir / f"{basename}_hillshade_multi.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        dy, dx = np.gradient(dem)

        azimuts = [315, 45, 225, 135]
        altitude = 30
        hillshades = []

        for az in azimuts:
            az_rad = np.radians(az)
            alt_rad = np.radians(altitude)
            slope = np.arctan(np.sqrt(dx**2 + dy**2))
            aspect = np.arctan2(dy, dx)
            hs = np.sin(alt_rad) * np.sin(slope) + \
                 np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect)
            hillshades.append(np.clip(hs, 0, 1))

        combined = np.mean(hillshades, axis=0)
        _save_tif(output, combined, transform, crs)
        logger.info(f"  ✓ Hillshade terminé ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur hillshade: {e}", exc_info=True)
        return None


def generate_slope(dem_file, basename, vis_dir, resolution):
    """Generate slope map (degrees)."""
    logger.info("  → Pente (Slope)...")
    t0 = time.time()
    output = vis_dir / f"{basename}_slope.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        dy, dx = np.gradient(dem)
        slope = np.arctan(np.sqrt(dx**2 + dy**2)) * 180 / np.pi
        _save_tif(output, slope, transform, crs)
        logger.info(f"  ✓ Pente terminée ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur slope: {e}", exc_info=True)
        return None


def generate_aspect(dem_file, basename, vis_dir, resolution):
    """Generate aspect (slope orientation) map."""
    logger.info("  → Aspect (Orientation)...")
    t0 = time.time()
    output = vis_dir / f"{basename}_aspect.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        dy, dx = np.gradient(dem)
        aspect = np.arctan2(dy, dx) * 180 / np.pi
        aspect = np.mod(aspect, 360)
        _save_tif(output, aspect, transform, crs)
        logger.info(f"  ✓ Aspect terminé ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur aspect: {e}", exc_info=True)
        return None


def generate_curvature(dem_file, basename, vis_dir, resolution):
    """Generate curvature (terrain concavity/convexity) map."""
    logger.info("  → Courbure (Curvature)...")
    t0 = time.time()
    output = vis_dir / f"{basename}_curvature.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        dz_dx = np.gradient(dem, axis=1)
        dz_dy = np.gradient(dem, axis=0)
        d2z_dx2 = np.gradient(dz_dx, axis=1)
        d2z_dy2 = np.gradient(dz_dy, axis=0)
        curvature = (d2z_dx2 + d2z_dy2) / 2
        _save_tif(output, curvature, transform, crs)
        logger.info(f"  ✓ Courbure terminée ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur curvature: {e}", exc_info=True)
        return None


def generate_solar(dem_file, basename, vis_dir, resolution):
    """Generate solar irradiance simulation."""
    logger.info("  → Éclairage Solaire...")
    t0 = time.time()
    output = vis_dir / f"{basename}_solar.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        dy, dx = np.gradient(dem)
        slope = np.arctan(np.sqrt(dx**2 + dy**2))
        aspect = np.arctan2(dy, dx)
        az_rad = np.radians(90)
        alt_rad = np.radians(30)
        solar = np.sin(alt_rad) * np.sin(slope) + \
                np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect)
        solar = np.clip(solar, 0, 1)
        _save_tif(output, solar, transform, crs)
        logger.info(f"  ✓ Solaire terminé ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur solar: {e}", exc_info=True)
        return None


# ============================================================
# GPU-accelerated visualizations
# ============================================================

def generate_lrm(dem_file, basename, vis_dir, resolution):
    """Local Relief Model - deviation from local mean (GPU if available)."""
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → Local Relief Model{gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_lrm.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        dem = to_gpu(dem_np)
        local_mean = xp_gaussian_filter(dem, sigma=15/resolution)
        lrm = dem - local_mean
        lrm_np = to_cpu(lrm).astype(np.float32)
        _save_tif(output, lrm_np, transform, crs)
        logger.info(f"  ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur LRM: {e}", exc_info=True)
        return None


def generate_svf(dem_file, basename, vis_dir, resolution):
    """Sky-View Factor - ray-tracing on 16 azimuths (GPU if available).

    For each pixel, trace rays in N directions, find the max horizon
    angle in each direction, then SVF = (1/N) * sum(cos²(horizon_angle)).
    Valleys/crevices have low SVF (obstructed sky), ridges/peaks have high SVF.
    """
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → Sky-View Factor (ray-tracing){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_svf.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        rows, cols = dem_np.shape
        res = resolution

        dem = to_gpu(dem_np)

        n_dirs = 16
        angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
        dx = np.cos(angles)
        dy = np.sin(angles)
        max_dist = int(50 / res)

        padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
        svf = xp.zeros_like(dem)

        for d_idx in range(n_dirs):
            ddx, ddy = dx[d_idx], dy[d_idx]
            horizon = xp.zeros_like(dem)

            for step in range(1, max_dist + 1):
                px = int(round(ddx * step))
                py = int(round(ddy * step))
                dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
                if dist_m < res * 0.5:
                    continue

                elev_diff = padded[max_dist + py:max_dist + py + rows,
                                  max_dist + px:max_dist + px + cols] - dem
                angle = xp.arctan2(elev_diff, dist_m)
                horizon = xp.where(xp.isnan(angle), horizon,
                                  xp.maximum(horizon, xp.nan_to_num(angle, nan=0)))

            svf += xp.cos(xp.pi / 2 - horizon) ** 2

        svf /= n_dirs
        svf_np = to_cpu(svf).astype(np.float32)
        _save_tif(output, svf_np, transform, crs)
        logger.info(f"  ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur SVF: {e}", exc_info=True)
        return None


def generate_openness(dem_file, basename, vis_dir, resolution, positive=True):
    """Positive/Negative Openness - true zenith/nadir angle computation (GPU if available).

    For each pixel, in 8 directions (N, NE, E, SE, S, SW, W, NW):
    - Positive openness: max zenith angle (angle from vertical to highest visible terrain)
    - Negative openness: max nadir angle (angle from vertical down to lowest terrain)
    Result is averaged across all 8 directions.
    """
    name = "positive_openness" if positive else "negative_openness"
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → {name.replace('_', ' ').title()} (ray-tracing){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_{name}.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        rows, cols = dem_np.shape
        res = resolution

        dem = to_gpu(dem_np)

        n_dirs = 8
        angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
        dx = np.cos(angles)
        dy = np.sin(angles)
        max_dist = int(50 / res)

        padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
        openness_sum = xp.zeros_like(dem)

        for d_idx in range(n_dirs):
            ddx, ddy = dx[d_idx], dy[d_idx]
            max_angle = xp.zeros_like(dem)

            for step in range(1, max_dist + 1):
                px = int(round(ddx * step))
                py = int(round(ddy * step))
                dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
                if dist_m < res * 0.5:
                    continue

                elev_diff = padded[max_dist + py:max_dist + py + rows,
                                  max_dist + px:max_dist + px + cols] - dem

                if positive:
                    angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m)
                else:
                    angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m)

                max_angle = xp.where(xp.isnan(angle), max_angle,
                                    xp.maximum(max_angle, xp.nan_to_num(angle, nan=0)))

            openness_sum += max_angle

        openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32)
        _save_tif(output, openness_result, transform, crs)
        logger.info(f"  ✓ {name} terminé ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur openness: {e}", exc_info=True)
        return None


def generate_mslrm(dem_file, basename, vis_dir, resolution):
    """Multi-Scale Relief Model (MSRM) - LRM at 5 scales combined (GPU if available)."""
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → Multi-Scale Relief Model (MSRM){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_mslrm.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        dem = to_gpu(dem_np)

        sigmas = [5, 10, 25, 50, 100]
        lrm_stack = []

        for sigma in sigmas:
            sigma_px = sigma / resolution
            local_mean = xp_gaussian_filter(dem, sigma=sigma_px)
            lrm = dem - local_mean
            lrm_norm = lrm / max(float(xp.nanstd(lrm)), 0.01)
            lrm_stack.append(lrm_norm)

        mslrm = xp.sqrt(xp.mean(xp.array(lrm_stack) ** 2, axis=0))
        mslrm_np = to_cpu(mslrm).astype(np.float32)
        _save_tif(output, mslrm_np, transform, crs)
        logger.info(f"  ✓ MSRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur MSRM: {e}", exc_info=True)
        return None


def generate_tpi(dem_file, basename, vis_dir, resolution):
    """Multi-Scale Topographic Position Index (GPU if available).

    TPI = elevation - mean(neighborhood).
    Computed at fine (5m) and broad (100m) scales.
    """
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → TPI multi-échelle{gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_tpi.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        dem = to_gpu(dem_np)

        fine_size = int(5 / resolution)
        if fine_size % 2 == 0:
            fine_size += 1
        tpi_fine = dem - xp_uniform_filter(dem, size=fine_size)

        broad_size = int(100 / resolution)
        if broad_size % 2 == 0:
            broad_size += 1
        tpi_broad = dem - xp_uniform_filter(dem, size=broad_size)

        fine_std = max(float(xp.nanstd(tpi_fine)), 0.01)
        broad_std = max(float(xp.nanstd(tpi_broad)), 0.01)
        tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std)
        tpi_np = to_cpu(tpi_combined).astype(np.float32)

        _save_tif(output, tpi_np, transform, crs)
        logger.info(f"  ✓ TPI terminé ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur TPI: {e}", exc_info=True)
        return None


# ============================================================
# Depression / hydrology
# ============================================================

def generate_depressions(dem_file, basename, vis_dir, resolution):
    """Depression detection using hydrological sink filling."""
    logger.info("  → Détection dépressions (hydrologique)...")
    t0 = time.time()
    output = vis_dir / f"{basename}_depressions.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)

        from scipy.ndimage import binary_dilation, generate_binary_structure

        dem_filled = dem.copy()
        nodata_mask = np.isnan(dem_filled)
        dem_filled[nodata_mask] = np.nanmax(dem) + 1000

        struct = generate_binary_structure(2, 2)
        changed = True
        iterations = 0
        max_iter = 100

        while changed and iterations < max_iter:
            from scipy.ndimage import minimum_filter as scipy_min_filter
            neighbor_min = scipy_min_filter(dem_filled, footprint=struct)
            sinks = (dem_filled < neighbor_min) & ~nodata_mask

            if not np.any(sinks):
                break

            new_dem = np.maximum(dem_filled, neighbor_min)
            new_dem[nodata_mask] = np.nan
            changed = np.any(new_dem != dem_filled)
            dem_filled = new_dem
            iterations += 1

        depressions = dem_filled - dem
        depressions[nodata_mask] = np.nan
        depressions = np.where(depressions > 0.01, depressions, 0)

        _save_tif(output, depressions, transform, crs)
        logger.info(f"  ✓ Dépressions terminé ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur dépressions: {e}", exc_info=True)
        return None


# ============================================================
# SAILORE
# ============================================================

def generate_sailore(dem_file, basename, vis_dir, resolution):
    """SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available).

    Kernel size adapts to local slope: flat areas get larger kernels,
    steep areas get smaller kernels.
    """
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → SAILORE (LRM adaptatif){gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_sailore.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        dem = to_gpu(dem_np)

        gy, gx = xp.gradient(dem, resolution)
        slope = xp.arctan(xp.sqrt(gx**2 + gy**2))
        slope_deg = xp.degrees(slope)

        sigma_min = 2.0 / resolution
        sigma_max = 25.0 / resolution
        slope_norm = xp.clip(slope_deg / 30.0, 0, 1)
        adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min)

        lrm_fine = dem - xp_gaussian_filter(dem, sigma=sigma_min)
        lrm_medium = dem - xp_gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2)
        lrm_coarse = dem - xp_gaussian_filter(dem, sigma=sigma_max)

        w_fine = slope_norm
        w_medium = 1 - 2 * xp.abs(slope_norm - 0.5)
        w_coarse = 1 - slope_norm
        w_total = w_fine + w_medium + w_coarse
        w_total[w_total == 0] = 1

        sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total
        sailore_np = to_cpu(sailore).astype(np.float32)

        _save_tif(output, sailore_np, transform, crs)
        logger.info(f"  ✓ SAILORE terminé ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur SAILORE: {e}", exc_info=True)
        return None


# ============================================================
# Roughness
# ============================================================

def generate_roughness(dem_file, basename, vis_dir, resolution):
    """Surface roughness - standard deviation of elevation in a window."""
    logger.info("  → Rugosité de surface...")
    t0 = time.time()
    output = vis_dir / f"{basename}_roughness.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        window_size = int(5 / resolution)
        if window_size % 2 == 0:
            window_size += 1

        def std_filter(arr):
            return np.nanstd(arr)

        roughness = generic_filter(dem.astype(np.float64), std_filter,
                                   size=window_size, mode='nearest')
        _save_tif(output, roughness, transform, crs)
        logger.info(f"  ✓ Rugosité terminée ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur rugosité: {e}", exc_info=True)
        return None


# ============================================================
# Anomalies
# ============================================================

def generate_anomalies(dem_file, basename, vis_dir, resolution):
    """Statistical anomaly detection - z-score of local relief + Local Moran's I."""
    logger.info("  → Détection anomalies statistiques...")
    t0 = time.time()
    output = vis_dir / f"{basename}_anomalies.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)

        lrm = dem - gaussian_filter(dem, sigma=15 / resolution)
        lrm_mean = np.nanmean(lrm)
        lrm_std = max(np.nanstd(lrm), 0.01)
        z_score = (lrm - lrm_mean) / lrm_std

        window = int(10 / resolution)
        if window % 2 == 0:
            window += 1

        local_mean = uniform_filter(z_score, size=window)
        morans_i = z_score * (local_mean - np.nanmean(z_score)) / np.nanstd(z_score)
        anomaly_score = np.abs(z_score) * np.sign(morans_i)

        _save_tif(output, anomaly_score, transform, crs)
        logger.info(f"  ✓ Anomalies terminé ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur anomalies: {e}", exc_info=True)
        return None


# ============================================================
# Wavelet
# ============================================================

def generate_wavelet(dem_file, basename, vis_dir, resolution):
    """Mexican Hat wavelet multi-scale analysis (GPU if available).

    CWT 2D at multiple scales to detect circular features.
    """
    gpu_tag = " [GPU]" if HAS_GPU else ""
    logger.info(f"  → Ondelette Mexican Hat multi-échelle{gpu_tag}...")
    t0 = time.time()
    output = vis_dir / f"{basename}_wavelet.tif"

    try:
        dem_np, transform, crs = _read_dem(dem_file)
        dem = to_gpu(dem_np)

        scales = [2, 5, 10, 20, 50]
        wavelet_stack = []

        for scale_m in scales:
            sigma_px = scale_m / resolution
            if HAS_GPU:
                from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace
                response = -gpu_gaussian_laplace(dem, sigma=sigma_px)
            else:
                from scipy.ndimage import gaussian_laplace
                response = to_gpu(-gaussian_laplace(to_cpu(dem).astype(np.float64), sigma=sigma_px))
            response /= max(float(xp.nanstd(response)), 0.01)
            wavelet_stack.append(response)

        combined = xp.sqrt(xp.mean(xp.array(wavelet_stack) ** 2, axis=0))
        combined_np = to_cpu(combined).astype(np.float32)

        _save_tif(output, combined_np, transform, crs)
        logger.info(f"  ✓ Ondelette terminée ({time.time()-t0:.1f}s){gpu_tag}")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur ondelette: {e}", exc_info=True)
        return None


# ============================================================
# Texture GLCM
# ============================================================

def generate_texture(dem_file, basename, vis_dir, resolution):
    """GLCM texture analysis on hillshade - contrast, entropy, homogeneity."""
    logger.info("  → Texture GLCM...")
    t0 = time.time()
    output = vis_dir / f"{basename}_texture.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)

        gy, gx = np.gradient(dem, resolution)
        slope = np.arctan(np.sqrt(gx**2 + gy**2))
        alt_rad = np.radians(45)
        az_rad = np.radians(315)
        shading = (np.sin(alt_rad) * np.cos(slope) +
                   np.cos(alt_rad) * np.sin(slope) *
                   np.cos(az_rad - np.arctan2(gy, gx)))
        hillshade = np.clip(shading, 0, 1)

        valid = hillshade[~np.isnan(hillshade)]
        if len(valid) == 0:
            raise ValueError("No valid data for texture analysis")
        lo, hi = np.percentile(valid, (1, 99))
        img = np.clip((hillshade - lo) / max(hi - lo, 0.001), 0, 1)
        img_uint8 = (img * 255).astype(np.uint8)

        window = int(5 / resolution)
        if window % 2 == 0:
            window += 1

        def local_variance(arr):
            return np.var(arr.astype(np.float64))

        def local_entropy(arr):
            hist, _ = np.histogram(arr.astype(np.float64), bins=16, range=(0, 256))
            hist = hist / max(hist.sum(), 1)
            hist = hist[hist > 0]
            return -np.sum(hist * np.log2(hist))

        def local_homogeneity(arr):
            arr_f = arr.astype(np.float64)
            return np.mean(1.0 / (1.0 + (arr_f - np.mean(arr_f)) ** 2))

        contrast = generic_filter(img_uint8.astype(np.float64), local_variance,
                                  size=window, mode='nearest')
        entropy = generic_filter(img_uint8.astype(np.float64), local_entropy,
                                 size=window, mode='nearest')
        homogeneity = generic_filter(img_uint8.astype(np.float64), local_homogeneity,
                                     size=window, mode='nearest')

        def norm(arr):
            valid_arr = arr[~np.isnan(arr)]
            if len(valid_arr) == 0:
                return arr
            std_val = max(np.std(valid_arr), 0.01)
            return (arr - np.mean(valid_arr)) / std_val

        texture_combined = (0.4 * norm(contrast) +
                           0.4 * norm(entropy) -
                           0.2 * norm(homogeneity))

        _save_tif(output, texture_combined, transform, crs)
        logger.info(f"  ✓ Texture terminée ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur texture GLCM: {e}", exc_info=True)
        return None


# ============================================================
# Flow accumulation
# ============================================================

def generate_flow(dem_file, basename, vis_dir, resolution):
    """Flow accumulation using D8 algorithm.

    Identifies drainage patterns, ditches, and enclosure boundaries.
    """
    logger.info("  → Accumulation de flux D8...")
    t0 = time.time()
    output = vis_dir / f"{basename}_flow.tif"

    try:
        dem, transform, crs = _read_dem(dem_file)
        rows, cols = dem.shape
        nodata_mask = np.isnan(dem)

        from scipy.ndimage import minimum_filter as scipy_min_filter, generate_binary_structure

        dem_filled = dem.copy()
        dem_filled[nodata_mask] = np.nanmax(dem) + 1000

        struct = generate_binary_structure(2, 2)
        for _ in range(50):
            neighbor_min = scipy_min_filter(dem_filled, footprint=struct)
            sinks = (dem_filled < neighbor_min) & ~nodata_mask
            if not np.any(sinks):
                break
            dem_filled = np.where(sinks, neighbor_min, dem_filled)

        dem_filled[nodata_mask] = np.nan

        dx8 = [1, 1, 0, -1, -1, -1, 0, 1]
        dy8 = [0, 1, 1, 1, 0, -1, -1, -1]
        dist8 = [1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2)]

        flow_dir = np.full((rows, cols), -1, dtype=np.int8)
        max_slope = np.full((rows, cols), 0.0)

        padded = np.pad(dem_filled, 1, mode='constant',
                       constant_values=np.nanmax(dem_filled) + 10000)

        for d in range(8):
            nx = 1 + dx8[d]
            ny = 1 + dy8[d]
            neighbor_elev = padded[ny:ny + rows, nx:nx + cols]
            slope = (dem_filled - neighbor_elev) / (dist8[d] * resolution)
            slope[nodata_mask] = -1
            better = slope > max_slope
            flow_dir[better] = d
            max_slope[better] = slope[better]

        flat_dem = dem_filled[~nodata_mask].flatten()
        valid_indices = np.where(~nodata_mask.flatten())[0]
        sort_order = valid_indices[np.argsort(-flat_dem)]

        flow_acc = np.ones((rows, cols), dtype=np.float32)
        flow_acc[nodata_mask] = 0

        for idx in sort_order:
            r, c = divmod(idx, cols)
            d = flow_dir[r, c]
            if d < 0:
                continue
            nr, nc = r + dy8[d], c + dx8[d]
            if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]:
                flow_acc[nr, nc] += flow_acc[r, c]

        flow_log = np.log1p(flow_acc)
        _save_tif(output, flow_log, transform, crs)
        logger.info(f"  ✓ Flux terminé ({time.time()-t0:.1f}s)")
        return output
    except Exception as e:
        logger.error(f"  ✗ Erreur flux: {e}", exc_info=True)
        return None