Improve visualizations: adaptive scales, revert z-score to std normalization

- MSRM/TPI/roughness/anomalies: revert z-score (x-mean)/std to std normalization x/std
  to preserve contrast and visibility of linear features (paths, ditches, trenches)
- MSRM: adaptive scales based on resolution, archaeological weight combination
- TPI: extend from 2 to 4 scales (3m/15m/50m/200m) with weighted combination
- Hillshade: 8 directions instead of 4, altitude 35° instead of 30°
- LRM: adaptive sigma based on resolution
- Openness: doubled radius (100m instead of 50m)
- Roughness: multi-scale (3m fine + 15m broad) instead of single 5x5 window
- Anomalies: uses MSRM multi-scale relief instead of single LRM 15m
- Wavelet: 8 adaptive scales, std normalization, archaeological weights
- Remove svf (Sky-View Factor) and local_dominance visualizations
- Add AVIF format support (default), quality 98
- Add multi-resolution support (-r 0.5,0.2)
- Improve Ctrl+C handling for immediate process termination
- Update rendering.py descriptions for all modified visualizations

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

CLAUDE.md

@@ -4,7 +4,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 
 ## Project Overview
 
-LiDAR archaeological processing pipeline that generates 20 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
+LiDAR archaeological processing pipeline that generates 16 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
 
 ## Commands
 
@@ -14,10 +14,15 @@ All commands run inside Docker. Use `./run.sh` as the primary interface.
 ./run.sh -g                                    # Standard run with GPU
 ./run.sh -g -w 4                               # GPU + 4 parallel workers
 ./run.sh -g -r 0.2                             # High resolution (0.2m/px)
+./run.sh -g -r 0.5,0.2                         # Multi-resolution (0.5m + 0.2m)
 ./run.sh --test                                 # Run unit tests
 ./run.sh -g --file LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc  # Single file
 ./run.sh --ground-classification csf            # Force CSF ground classification (complex terrain)
 ./run.sh -g --keep-tif                          # Keep TIFF files (allows WebP regeneration without recalculating DTM)
+./run.sh -g --only hillshade svf lrm            # Only generate specific visualizations
+./run.sh -g --skip ortho topo                   # Exclude specific visualizations
+./run.sh -g --quality 90                       # WebP quality 90 (default: 85)
+./run.sh -g --lossless                          # Lossless WebP compression
 ./run.sh                                        # Print help (no args)
 ```
 
@@ -32,12 +37,12 @@ docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data
 ### Module responsibilities
 
 - **`cli.py`** — argparse + logging setup. Entry point via `python -m lidar_pipeline`.
-- **`pipeline.py`** — `LidarArchaeoPipeline` orchestrator. `VIZ_STEPS` registry maps names to generate functions. `FilePrefixFilter` for parallel logging. Creates `SharedDEM` once per file and passes it to all visualizations.
-- **`dtm.py`** — PDAL ground classification (SMRF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`.
-- **`visualizations.py`** — 15 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution, shared=None)` and return a TIF path or None. `SharedDEM` class pre-computes gradient, NaN mask, LRM to avoid redundant I/O and computation.
+- **`pipeline.py`** — `LidarArchaeoPipeline` orchestrator. `VIZ_STEPS` registry maps names to generate functions. `FilePrefixFilter` for parallel logging. Creates `SharedDEM` once per file and passes it to all visualizations. Multi-resolution support: `self.resolutions` list, `_res_suffix()` for naming, `generate_all_visualizations()` accepts `vis_dir` override.
+- **`dtm.py`** — PDAL ground classification (SMRF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`. `create_dtm_fast()` accepts `output_suffix` for multi-resolution DTM naming.
+- **`visualizations.py`** — 13 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution, shared=None)` and return a TIF path or None. `SharedDEM` class pre-computes gradient, NaN mask, LRM to avoid redundant I/O and computation. Lazy evaluation: properties computed on first access.
 - **`gpu.py`** — CuPy/numpy abstraction: `HAS_GPU`, `to_gpu()`, `to_cpu()`, `xp_gaussian_filter()`, `xp_uniform_filter()`, `xp_minimum_filter()`, `gpu_cleanup()`. Falls back to CPU gracefully.
 - **`ign.py`** — IGN WMTS tile download + overlay generation for orthophoto and topographic maps.
-- **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. `generate_pdf_report()` creates A3 PDF.
+- **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. Quality parameter controls WebP compression (default 85).
 
 ### SharedDEM optimization
 
@@ -73,6 +78,10 @@ DTM small gaps (< 1m from existing data) are filled using `rasterio.fill.fillnod
 
 Uses priority-flood algorithm (Wang & Liu 2006) for sink filling, which is O(n log n) instead of iterative minimum_filter. D8 accumulation uses numba JIT; falls back to pure Python if numba unavailable.
 
+### Multi-resolution
+
+`-r 0.5,0.2` processes each tile at both 0.5m and 0.2m. Ground classification is shared (done once per tile). Each resolution gets its own DTM (`_dtm.tif` / `_dtm_r0p2.tif`) and visualization subdirectory (`basename/` / `basename_r0p2/`).
+
 ### Parallel processing
 
 Uses `ProcessPoolExecutor` with `'spawn'` start method (required for CUDA). Each worker gets its own temp directory (`temp_{basename}`). `_process_file_standalone()` configures its own logger with `_file_filter` for per-file log prefixes.
@@ -82,6 +91,6 @@ Uses `ProcessPoolExecutor` with `'spawn'` start method (required for CUDA). Each
 - **Language**: UI messages and comments in French. Code identifiers in English.
 - **Logging**: Use `logger = logging.getLogger("lidar")`. Prefix per-file logs via `_file_filter.basename`.
 - **GPU pattern**: `arr_gpu = to_gpu(arr)` → compute → `result = to_cpu(arr_gpu)` → `gpu_cleanup()` between visualizations.
-- **Output format**: Visualizations saved as WebP. TIFF intermediates deleted by default. Use `--keep-tif` to keep DTM+TIF for WebP regeneration with `--force`. No COGs or viewer.
+- **Output format**: Visualizations saved as AVIF (quality 98 by default, best quality/size ratio). Use `--format webp` for WebP output. TIFF intermediates deleted by default. Use `--keep-tif` to keep DTM+TIF for regeneration with `--force`. No PDF reports, no COGs or viewer.
 - **Compression**: TIF intermediates use `deflate` compression (faster than LZW for float32 data).
 - **Tests**: Run only inside Docker via `./run.sh --test`. Synthetic DEM fixture in `tests/conftest.py`.

Dockerfile

@@ -41,7 +41,8 @@ RUN pip3 install --no-cache-dir \
     titiler.core \
     fastapi \
     uvicorn \
-    piexif
+    piexif \
+    pillow-avif-plugin
 
 # Install CuPy for GPU acceleration (optional - will fallback to numpy if not available)
 RUN pip3 install --no-cache-dir cupy-cuda12x || echo "CuPy not available - GPU acceleration disabled"

download_lidar.sh

@@ -0,0 +1,584 @@
+#!/bin/bash
+# download_lidar.sh — Téléchargement de fichiers LiDAR HD depuis l'IGN
+#
+# Usage:
+#   ./download_lidar.sh LHD_FXX_1049_6895_PTS_LAMB93_IGN69.copc.laz
+#   ./download_lidar.sh 1049_6895
+#   ./download_lidar.sh 1049_6895 1029_6884
+#   ./download_lidar.sh --list-zones
+#   ./download_lidar.sh --search 1049
+#   ./download_lidar.sh --zone RE --output input/
+#   ./download_lidar.sh --zone RE --dry-run
+#
+# L'API IGN expose les fichiers organisés par zones (RE, SE, AE, etc.)
+# Le script cherche automatiquement dans quelle zone se trouve chaque tuile.
+
+set -euo pipefail
+
+BASE_URL="https://data.geopf.fr/telechargement/download/LiDARHD-NUALID"
+API_URL="https://data.geopf.fr/telechargement/resource/LiDARHD-NUALID"
+OUTPUT_DIR="${OUTPUT_DIR:-.}"
+
+# Colors
+RED='\033[0;31m'
+GREEN='\033[0;32m'
+YELLOW='\033[1;33m'
+CYAN='\033[0;36m'
+NC='\033[0m'
+
+log_info()  { echo -e "${CYAN}[INFO]${NC} $*"; }
+log_ok()    { echo -e "${GREEN}[OK]${NC} $*"; }
+log_warn()  { echo -e "${YELLOW}[WARN]${NC} $*"; }
+log_err()   { echo -e "${RED}[ERR]${NC} $*"; }
+
+# Normalize a tile identifier to a full filename
+# Supports: "1049_6895", "LHD_FXX_1049_6895_PTS_LAMB93_IGN69", "LHD_FXX_1049_6895_PTS_LAMB93_IGN69.copc.laz"
+normalize_name() {
+    local input="$1"
+    # If it's already a full filename, return as-is
+    if [[ "$input" == LHD_FXX_*_PTS_LAMB93_IGN69.copc.laz ]]; then
+        echo "$input"
+        return
+    fi
+    if [[ "$input" == LHD_FXX_*_PTS_LAMB93_IGN69 ]]; then
+        echo "${input}.copc.laz"
+        return
+    fi
+    # Bare coordinates like "1049_6895"
+    echo "LHD_FXX_${input}_PTS_LAMB93_IGN69.copc.laz"
+}
+
+# List all available zone releases
+list_zones() {
+    log_info "Récupération de la liste des zones..."
+    local page=1
+    local total_pages=""
+
+    while true; do
+        local content
+        content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
+
+        if [ -z "$content" ]; then
+            break
+        fi
+
+        # Parse zone codes and dates
+        echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    for entry in root.findall('atom:entry', ns):
+        title = entry.find('atom:title', ns)
+        if title is not None and title.text:
+            # Extract zone code and date
+            parts = title.text.split('_')
+            # NUALHD_1-0__LAZ_LAMB93_{ZONE}_{DATE}
+            zone = parts[-2] if len(parts) >= 2 else '?'
+            date = parts[-1] if len(parts) >= 1 else '?'
+            print(f'{zone}\t{date}\t{title.text}')
+except Exception:
+    pass
+" 2>/dev/null || true
+
+        # Check if more pages
+        if [ -z "$total_pages" ]; then
+            total_pages=$(echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom', 'gpf_dl': 'https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    pc = root.get('{https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd}pagecount', '1')
+    print(pc)
+except:
+    print('1')
+" 2>/dev/null)
+        fi
+
+        if [ "$page" -ge "${total_pages:-1}" ]; then
+            break
+        fi
+        page=$((page + 1))
+    done | sort
+}
+
+# Search for a tile across all zones
+search_tile() {
+    local tile_name="$1"
+    local filename
+    filename=$(normalize_name "$tile_name")
+
+    log_info "Recherche de ${filename} dans toutes les zones..."
+
+    local page=1
+    local found=0
+    local total_pages=""
+
+    while true; do
+        local content
+        content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
+        [ -z "$content" ] && break
+
+        # Get sub-resources and search for the file
+        echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    for entry in root.findall('atom:entry', ns):
+        title = entry.find('atom:title', ns)
+        if title is not None and title.text:
+            print(title.text)
+except:
+    pass
+" 2>/dev/null | while read -r release; do
+            # Check if file exists in this release
+            local url="${BASE_URL}/${release}/${filename}"
+            local code
+            code=$(curl -sI -L -o /dev/null -w "%{http_code}" "$url" 2>/dev/null)
+            if [ "$code" = "200" ]; then
+                log_ok "Trouvé: ${filename} dans la zone ${release}"
+                echo "  URL: ${url}"
+                found=1
+            fi
+        done
+
+        if [ -z "$total_pages" ]; then
+            total_pages=$(echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom', 'gpf_dl': 'https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    pc = root.get('{https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd}pagecount', '1')
+    print(pc)
+except:
+    print('1')
+" 2>/dev/null)
+        fi
+
+        if [ "$page" -ge "${total_pages:-1}" ]; then
+            break
+        fi
+        page=$((page + 1))
+    done
+
+    if [ "$found" -eq 0 ]; then
+        log_err "${filename} non trouvé dans aucune zone"
+    fi
+}
+
+# Download a single file with resume support
+download_file() {
+    local url="$1"
+    local output_path="$2"
+    local filename
+    filename=$(basename "$output_path")
+
+    # Resume download with retries
+    local max_retries=15
+    local retry=0
+
+    while [ $retry -lt $max_retries ]; do
+        retry=$((retry + 1))
+        local result
+        result=$(curl -L --http1.1 -C - -o "$output_path" "$url" 2>&1)
+        local exit_code=$?
+
+        if [ $exit_code -eq 0 ]; then
+            # Verify it's a valid LAS/LAZ file
+            local magic
+            magic=$(head -c 4 "$output_path" 2>/dev/null | xxd -p 2>/dev/null || echo "")
+            if [ "$magic" = "4c415346" ]; then
+                local size
+                size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
+                log_ok "${filename} téléchargé (${size} octets)"
+                return 0
+            else
+                log_err "${filename} n'est pas un fichier LAZ valide (magic: ${magic})"
+                rm -f "$output_path"
+                return 1
+            fi
+        fi
+
+        # curl exit code 18 = partial download, resume
+        if [ $exit_code -eq 18 ]; then
+            local size
+            size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
+            log_warn "Téléchargement partiel (${size} octets), tentative ${retry}/${max_retries}..."
+            sleep 2
+            continue
+        fi
+
+        # HTTP 404
+        if echo "$result" | grep -q "404"; then
+            log_err "Fichier non trouvé (404): ${url}"
+            rm -f "$output_path"
+            return 1
+        fi
+
+        log_err "Erreur curl (${exit_code}): ${result}"
+        return 1
+    done
+
+    # If we exhausted retries, check if partial file is usable
+    local magic
+    magic=$(head -c 4 "$output_path" 2>/dev/null | xxd -p 2>/dev/null || echo "")
+    if [ "$magic" = "4c415346" ]; then
+        local size
+        size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
+        log_warn "Fichier partiel mais valide: ${filename} (${size} octets)"
+        return 0
+    fi
+
+    log_err "Échec après ${max_retries} tentatives"
+    rm -f "$output_path"
+    return 1
+}
+
+# Find which zone contains a given tile and download it
+find_and_download() {
+    local tile_name="$1"
+    local output_dir="$2"
+    local filename
+    filename=$(normalize_name "$tile_name")
+    local output_path="${output_dir}/${filename}"
+
+    # If file already exists and is valid, skip
+    if [ -f "$output_path" ]; then
+        local magic
+        magic=$(head -c 4 "$output_path" 2>/dev/null | xxd -p 2>/dev/null || echo "")
+        if [ "$magic" = "4c415346" ]; then
+            local size
+            size=$(stat -c%s "$output_path" 2>/dev/null || echo "0")
+            if [ "$size" -gt 1000000 ]; then
+                log_info "${filename} déjà présent (${size} octets) — ignoré"
+                return 0
+            fi
+        fi
+        # File exists but seems invalid/truncated, re-download
+        log_warn "${filename} existe mais semble incomplet, re-téléchargement..."
+    fi
+
+    log_info "Recherche de ${filename}..."
+
+    # Try all zone releases
+    local page=1
+    local total_pages=""
+
+    while true; do
+        local content
+        content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
+        [ -z "$content" ] && break
+
+        echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    for entry in root.findall('atom:entry', ns):
+        title = entry.find('atom:title', ns)
+        if title is not None and title.text:
+            print(title.text)
+except:
+    pass
+" 2>/dev/null | while read -r release; do
+            local url="${BASE_URL}/${release}/${filename}"
+            local code
+            code=$(curl -sI -L -o /dev/null -w "%{http_code}" "$url" 2>/dev/null || echo "000")
+            if [ "$code" = "200" ]; then
+                log_ok "Trouvé dans la zone ${release}"
+                download_file "$url" "$output_path"
+                return 0
+            fi
+        done
+
+        if [ -z "$total_pages" ]; then
+            total_pages=$(echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom', 'gpf_dl': 'https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    pc = root.get('{https://data.geopf.fr/annexes/ressources/xsd/gpf_dl.xsd}pagecount', '1')
+    print(pc)
+except:
+    print('1')
+" 2>/dev/null)
+        fi
+
+        if [ "$page" -ge "${total_pages:-1}" ]; then
+            break
+        fi
+        page=$((page + 1))
+    done
+
+    log_err "${filename} non trouvé dans aucune zone disponible"
+    return 1
+}
+
+# Download all tiles from a specific zone
+download_zone() {
+    local zone="$1"
+    local output_dir="$2"
+    local dry_run="${3:-false}"
+
+    # Find the release name for this zone
+    local release=""
+    local page=1
+
+    while true; do
+        local content
+        content=$(curl -s "${API_URL}?page=${page}&limit=50") || true
+        [ -z "$content" ] && break
+
+        release=$(echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    for entry in root.findall('atom:entry', ns):
+        title = entry.find('atom:title', ns)
+        if title is not None and title.text:
+            # Extract zone code from title like NUALHD_1-0__LAZ_LAMB93_RE_2025-02-17
+            parts = title.text.split('_')
+            z = parts[-2] if len(parts) >= 2 else ''
+            if z == '${zone}':
+                print(title.text)
+                break
+except:
+    pass
+" 2>/dev/null)
+
+        if [ -n "$release" ]; then
+            break
+        fi
+        page=$((page + 1))
+        # Safety limit
+        if [ $page -gt 25 ]; then
+            break
+        fi
+    done
+
+    if [ -z "$release" ]; then
+        log_err "Zone '${zone}' non trouvée"
+        return 1
+    fi
+
+    log_info "Zone ${zone}: ${release}"
+
+    # Get list of files in this release
+    local sub_url="${API_URL}/${release}?limit=50"
+    local total_files=0
+    local page=1
+
+    while true; do
+        local content
+        content=$(curl -s "${API_URL}/${release}?page=${page}&limit=50") || true
+        [ -z "$content" ] && break
+
+        local files
+        files=$(echo "$content" | python3 -c "
+import sys, xml.etree.ElementTree as ET
+ns = {'atom': 'http://www.w3.org/2005/Atom'}
+try:
+    tree = ET.parse(sys.stdin)
+    root = tree.getroot()
+    for entry in root.findall('atom:entry', ns):
+        title = entry.find('atom:title', ns)
+        if title is not None and title.text:
+            print(title.text)
+except:
+    pass
+" 2>/dev/null)
+
+        if [ -z "$files" ]; then
+            break
+        fi
+
+        local count
+        count=$(echo "$files" | wc -l)
+        total_files=$((total_files + count))
+
+        if [ "$dry_run" = "true" ]; then
+            echo "$files" | head -10
+            echo "... (${count} fichiers sur cette page)"
+        else
+            echo "$files" | while read -r fname; do
+                local url="${BASE_URL}/${release}/${fname}"
+                download_file "$url" "${output_dir}/${fname}"
+            done
+        fi
+
+        page=$((page + 1))
+        # Safety limit
+        if [ $page -gt 500 ]; then
+            break
+        fi
+    done
+
+    log_info "Total: ${total_files} fichiers dans la zone ${zone}"
+}
+
+# Show usage
+usage() {
+    cat <<'EOF'
+download_lidar.sh — Téléchargement de fichiers LiDAR HD depuis l'IGN
+
+Usage:
+  ./download_lidar.sh <tuile> [tuile2...]    Télécharger une ou plusieurs tuiles
+  ./download_lidar.sh --search <motif>        Chercher une tuile dans toutes les zones
+  ./download_lidar.sh --list-zones            Lister les zones disponibles
+  ./download_lidar.sh --zone <ZONE> [--dry-run]  Télécharger une zone entière
+  ./download_lidar.sh --check [dossier/]      Vérifier l'intégrité des fichiers LAZ
+
+Formats de tuile acceptés:
+  1049_6895                                    Coordonnées seules
+  LHD_FXX_1049_6895_PTS_LAMB93_IGN69          Nom complet sans extension
+  LHD_FXX_1049_6895_PTS_LAMB93_IGN69.copc.laz Nom complet avec extension
+
+Options:
+  -o, --output DIR    Répertoire de destination (défaut: .)
+  -n, --dry-run       Afficher sans télécharger (avec --zone)
+  --check [DIR]       Vérifier l'intégrité des fichiers LAZ (défaut: ./)
+  --list-zones        Lister toutes les zones disponibles
+  --search MOTIF      Chercher un motif dans les noms de tuile
+  --zone ZONE         Télécharger tous les fichiers d'une zone
+  -h, --help          Afficher cette aide
+
+Exemples:
+  ./download_lidar.sh 1049_6895
+  ./download_lidar.sh LHD_FXX_0713_6347_PTS_LAMB93_IGN69.copc.laz -o input/
+  ./download_lidar.sh --search 1049
+  ./download_lidar.sh --list-zones
+  ./download_lidar.sh --zone RE -o input/
+  ./download_lidar.sh --check input/
+EOF
+}
+
+# Check integrity of LAZ files
+check_files() {
+    local dir="${1:-.}"
+    local ok=0
+    local fail=0
+
+    log_info "Vérification de l'intégrité des fichiers LAZ dans ${dir}/"
+
+    for f in "${dir}"/*.copc.laz "${dir}"/*.laz; do
+        [ -f "$f" ] || continue
+        local magic
+        magic=$(head -c 4 "$f" 2>/dev/null | xxd -p 2>/dev/null || echo "")
+        local size
+        size=$(stat -c%s "$f" 2>/dev/null || echo "0")
+
+        if [ "$magic" = "4c415346" ] && [ "$size" -gt 1000000 ]; then
+            log_ok "$(basename $f) (${size} octets)"
+            ok=$((ok + 1))
+        elif [ "$magic" = "7b227469" ]; then
+            log_err "$(basename $f) est une page HTML (404), pas un fichier LAZ"
+            fail=$((fail + 1))
+        elif [ "$size" -lt 1000 ]; then
+            log_err "$(basename $f) trop petit (${size} octets) — probablement corrompu"
+            fail=$((fail + 1))
+        else
+            log_warn "$(basename $f) magic=${magic}, taille=${size} — vérification nécessaire"
+            fail=$((fail + 1))
+        fi
+    done
+
+    log_info "Résultat: ${ok} OK, ${fail} problème(s)"
+    return 0
+}
+
+# Main
+CMD=""
+TILES=()
+ZONE=""
+DRY_RUN=false
+
+while [ $# -gt 0 ]; do
+    case "$1" in
+        -h|--help)
+            usage
+            exit 0
+            ;;
+        --list-zones)
+            CMD="list-zones"
+            shift
+            ;;
+        --search)
+            CMD="search"
+            shift
+            TILES+=("$1")
+            shift
+            ;;
+        --zone)
+            CMD="zone"
+            shift
+            ZONE="$1"
+            shift
+            ;;
+        --check)
+            CMD="check"
+            shift
+            OUTPUT_DIR="${1:-.}"
+            [ -d "$OUTPUT_DIR" ] && shift || true
+            ;;
+        -o|--output)
+            shift
+            OUTPUT_DIR="$1"
+            shift
+            ;;
+        -n|--dry-run)
+            DRY_RUN=true
+            shift
+            ;;
+        -*)
+            log_err "Option inconnue: $1"
+            usage
+            exit 1
+            ;;
+        *)
+            TILES+=("$1")
+            shift
+            ;;
+    esac
+done
+
+case "${CMD}" in
+    list-zones)
+        list_zones
+        ;;
+    search)
+        for tile in "${TILES[@]}"; do
+            search_tile "$tile"
+        done
+        ;;
+    zone)
+        if [ -z "$ZONE" ]; then
+            log_err "Zone non spécifiée. Utilisez --zone <CODE>"
+            exit 1
+        fi
+        mkdir -p "$OUTPUT_DIR"
+        download_zone "$ZONE" "$OUTPUT_DIR" "$DRY_RUN"
+        ;;
+    check)
+        check_files "$OUTPUT_DIR"
+        ;;
+    *)
+        # Default: download tiles
+        if [ ${#TILES[@]} -eq 0 ]; then
+            usage
+            exit 1
+        fi
+        mkdir -p "$OUTPUT_DIR"
+        for tile in "${TILES[@]}"; do
+            find_and_download "$tile" "$OUTPUT_DIR"
+        done
+        ;;
+esac

lidar_pipeline/cli.py

@@ -97,9 +97,9 @@ Exemples:
     )
     parser.add_argument(
         "-r", "--resolution",
-        type=float,
-        default=0.5,
-        help="Résolution en mètres par pixel (défaut: 0.5)"
+        type=str,
+        default="0.5",
+        help="Résolution en m/px, ou multiples séparées par virgules (défaut: 0.5, ex: 0.5,0.2)"
     )
     parser.add_argument(
         "-w", "--workers",
@@ -128,6 +128,37 @@ Exemples:
         default="auto",
         help="Méthode de classification du sol : auto (détection), smrf, csf (défaut: auto)"
     )
+    parser.add_argument(
+        "--quality",
+        type=int,
+        default=98,
+        help="Qualité image (1-100, défaut: 98). Utilisez 100 pour lossless."
+    )
+    parser.add_argument(
+        "--lossless",
+        action="store_true",
+        help="Forcer la compression lossless (équivalent à --quality 100)"
+    )
+    parser.add_argument(
+        "--format",
+        choices=["webp", "avif"],
+        default="avif",
+        help="Format de sortie : avif (défaut, meilleure qualité) ou webp"
+    )
+    parser.add_argument(
+        "--only",
+        nargs="+",
+        type=str,
+        default=None,
+        help="Générer uniquement ces visualisations (ex: --only hillshade svf lrm)"
+    )
+    parser.add_argument(
+        "--skip",
+        nargs="+",
+        type=str,
+        default=None,
+        help="Exclure ces visualisations (ex: --skip ortho topo)"
+    )
     parser.add_argument(
         "--file",
         nargs="+",
@@ -168,6 +199,14 @@ Exemples:
     log_gpu_status()
 
     try:
+        quality = 100 if args.lossless else args.quality
+        # Parse --only and --skip: accept comma-separated values
+        only_viz = None
+        if args.only:
+            only_viz = [v.strip() for item in args.only for v in item.split(',')]
+        skip_viz = None
+        if args.skip:
+            skip_viz = [v.strip() for item in args.skip for v in item.split(',')]
         pipeline = LidarArchaeoPipeline(
             input_dir=args.input,
             output_dir=args.output,
@@ -177,6 +216,10 @@ Exemples:
             ground_method=args.ground_classification,
             force_classify=args.force_classification,
             keep_tif=args.keep_tif,
+            quality=quality,
+            only_viz=only_viz,
+            skip_viz=skip_viz,
+            output_format=args.format,
         )
 
         # If --file is specified, process only matching files
@@ -187,10 +230,16 @@ Exemples:
             # Also supports bare name without .copc (e.g. LHD_FXX_1000_6882_PTS_LAMB93_IGN69)
             selected_files = []
             for pattern in args.file:
-                matches = (list(input_dir.glob(f"{pattern}.laz"))
-                          + list(input_dir.glob(f"{pattern}.las"))
-                          + list(input_dir.glob(f"{pattern}.copc.laz"))
-                          + list(input_dir.glob(f"{pattern}.copc.las")))
+                # Try exact filename first (e.g. LHD_FXX_...copc.laz)
+                exact_match = input_dir / pattern
+                if exact_match.exists() and exact_match.is_file():
+                    matches = [exact_match]
+                else:
+                    # Try with added extensions (e.g. pattern=LHD_FXX_...IGN69)
+                    matches = (list(input_dir.glob(f"{pattern}.laz"))
+                              + list(input_dir.glob(f"{pattern}.las"))
+                              + list(input_dir.glob(f"{pattern}.copc.laz"))
+                              + list(input_dir.glob(f"{pattern}.copc.las")))
                 # Remove duplicates
                 matches = list(dict.fromkeys(matches))
                 if not matches:
@@ -235,9 +284,15 @@ def _kill_orphan_pdal(signum=None, frame=None):
     """Kill orphan PDAL processes on interrupt or exit."""
     import subprocess
     try:
-        subprocess.run(["pkill", "-f", "pdal"], capture_output=True, timeout=5)
+        subprocess.run(["pkill", "-9", "-f", "pdal"], capture_output=True, timeout=3)
     except Exception:
         pass
     if signum is not None:
-        logger.info("Interruption — nettoyage des processus PDAL")
+        logger.info("Interruption — nettoyage des processus")
+        # Force-kill all child processes immediately
+        try:
+            import os
+            os.killpg(os.getpgrp(), signal.SIGKILL)
+        except Exception:
+            pass
         sys.exit(130)

lidar_pipeline/dtm.py

@@ -128,17 +128,21 @@ def validate_laz(laz_file):
     """Quick integrity check for a LAZ/LAS file.
 
     Tries laspy first (fast header read), then PDAL as fallback for COPC files
-    that laspy cannot read.
+    that laspy cannot read. Also checks that the file contains points.
 
     Returns:
-        True if file is readable, False otherwise.
+        True if file is readable and contains points, False otherwise.
     """
     # Try laspy first (fast)
     import laspy
     try:
         with laspy.open(str(laz_file)) as f:
             header = f.header
-            _ = header.point_count
+            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
         return True
     except Exception:
         pass
@@ -151,6 +155,17 @@ def validate_laz(laz_file):
             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
             return True
         logger.error(f"  ✗ Fichier illisible: {laz_file.name}")
         logger.error(f"    PDAL: {result.stderr.strip()[:200]}")
@@ -200,6 +215,9 @@ def _read_with_pdal(laz_file):
             os.unlink(tmp_path)
         except Exception:
             pass
+        if len(las.points) == 0:
+            logger.warning(f"  PDAL: conversion réussie mais 0 points")
+            return None
         return las
 
     except Exception as e:
@@ -238,6 +256,10 @@ def detect_ground_method(laz_file):
         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'
+
     z = np.array(las.z)
 
     # Height variance (always available)
@@ -321,6 +343,11 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
             ["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)
+            return None
         logger.info(f"  ✓ Classification sol {method.upper()} terminée")
         return output_las
     except subprocess.CalledProcessError as e:
@@ -375,7 +402,7 @@ def _repair_laz_with_laspy(input_laz, output_las):
         return False
 
 
-def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False):
+def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output_suffix=""):
     """Create DTM using fast binning method with gap filling.
 
     Args:
@@ -384,11 +411,12 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False):
         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).
 
     Returns:
         Path to output DTM GeoTIFF, or None on failure.
     """
-    output_tif = dtm_dir / f"{basename}_dtm.tif"
+    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}")
@@ -409,6 +437,10 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False):
             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])
@@ -451,7 +483,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False):
             logger.info(f"    {remaining:,} pixels restent sans données (grands écarts)")
 
         # Save as GeoTIFF
-        output_tif = dtm_dir / f"{basename}_dtm.tif"
+        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(

lidar_pipeline/pipeline.py

@@ -58,14 +58,14 @@ from .dtm import classify_ground, create_dtm_fast
 from .visualizations import (
     SharedDEM,
     generate_hillshade, generate_slope, generate_aspect, generate_curvature,
-    generate_lrm, generate_svf, generate_openness,
+    generate_lrm, generate_openness,
     generate_mslrm, generate_tpi, generate_sailore,
     generate_roughness, generate_anomalies, generate_wavelet,
-    generate_flow, generate_local_dominance,
+    generate_flow,
 )
 from .gpu import gpu_cleanup
 from .ign import generate_ign_overlay
-from .rendering import tif_to_png, generate_pdf_report
+from .rendering import tif_to_png
 
 
 # Ordered list of visualization steps.
@@ -76,7 +76,6 @@ VIZ_STEPS = [
     ('slope',      generate_slope),
     ('aspect',     generate_aspect),
     ('curvature',  generate_curvature),
-    ('svf',        generate_svf),
     ('lrm',        generate_lrm),
     ('pos_open',   lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=True, shared=shared)),
     ('neg_open',   lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=False, shared=shared)),
@@ -87,7 +86,6 @@ VIZ_STEPS = [
     ('anomalies',  generate_anomalies),
     ('wavelet',    generate_wavelet),
     ('flow',       generate_flow),
-    ('local_dominance', generate_local_dominance),
     ('ortho',      lambda d, b, v, r: generate_ign_overlay(
                     d, b, v, r,
                     layer='ORTHOIMAGERY.ORTHOPHOTOS',
@@ -108,15 +106,26 @@ VIZ_STEPS = [
 class LidarArchaeoPipeline:
     """Orchestrates the LiDAR archaeological analysis pipeline."""
 
-    def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False):
+    def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False, quality=98, only_viz=None, skip_viz=None, output_format='avif'):
         self.input_dir = Path(input_dir)
         self.output_dir = Path(output_dir)
-        self.resolution = resolution
+        # Accept single float or comma-separated string for multi-resolution
+        if isinstance(resolution, str):
+            self.resolutions = [float(r.strip()) for r in resolution.split(',')]
+        elif isinstance(resolution, (list, tuple)):
+            self.resolutions = [float(r) for r in resolution]
+        else:
+            self.resolutions = [float(resolution)]
+        self.resolution = self.resolutions[0]  # Primary resolution (backward compat)
         self.workers = workers
         self.force = force
         self.ground_method = ground_method
         self.force_classify = force_classify
         self.keep_tif = keep_tif
+        self.quality = quality
+        self.only_viz = only_viz
+        self.skip_viz = skip_viz
+        self.output_format = output_format
         self.temp_dir = self.output_dir / "temp"
 
         if not self.input_dir.exists():
@@ -127,20 +136,43 @@ class LidarArchaeoPipeline:
 
         self.dtm_dir = self.output_dir / "DTM"
         self.vis_dir = self.output_dir / "visualisations"
-        self.pdf_dir = self.output_dir / "rapports"
 
-        for d in [self.dtm_dir, self.vis_dir, self.pdf_dir]:
+        for d in [self.dtm_dir, self.vis_dir]:
             d.mkdir(exist_ok=True)
 
+        # Filter visualizations based on --only / --skip
+        all_viz_names = [name for name, _ in VIZ_STEPS]
+        if only_viz:
+            invalid = set(only_viz) - set(all_viz_names)
+            if invalid:
+                raise ValueError(f"Visualisations inconnues: {', '.join(invalid)}. Disponibles: {', '.join(all_viz_names)}")
+            self.viz_steps = [(n, f) for n, f in VIZ_STEPS if n in only_viz]
+        elif skip_viz:
+            invalid = set(skip_viz) - set(all_viz_names)
+            if invalid:
+                raise ValueError(f"Visualisations inconnues: {', '.join(invalid)}. Disponibles: {', '.join(all_viz_names)}")
+            self.viz_steps = [(n, f) for n, f in VIZ_STEPS if n not in skip_viz]
+        else:
+            self.viz_steps = VIZ_STEPS
+
         logger.info("Pipeline initialisé")
         logger.info(f"  Entrée      : {self.input_dir}")
         logger.info(f"  Sortie      : {self.output_dir}")
-        logger.info(f"  Résolution  : {resolution}m/px")
+        if len(self.resolutions) > 1:
+            logger.info(f"  Résolutions : {', '.join(f'{r}m/px' for r in self.resolutions)}")
+        else:
+            logger.info(f"  Résolution  : {self.resolution}m/px")
         logger.info(f"  Workers     : {workers}")
         logger.info(f"  Force       : {'OUI' if self.force else 'non (skip existing)'}")
         logger.info(f"  Classification sol : {self.ground_method}")
         logger.info(f"  Force classif.: {'OUI' if self.force_classify else 'non'}")
         logger.info(f"  Keep TIFF    : {'OUI' if self.keep_tif else 'non'}")
+        logger.info(f"  Qualité {self.output_format.upper()}: {self.quality if self.quality < 100 else 'lossless'}")
+        if only_viz:
+            logger.info(f"  Visualisations: uniquement {', '.join(only_viz)}")
+        elif skip_viz:
+            logger.info(f"  Visualisations: tout sauf {', '.join(skip_viz)}")
+        logger.info(f"  Visualisations: {len(self.viz_steps)}/{len(VIZ_STEPS)}")
 
     def find_laz_files(self):
         """Find all LAZ/LAS files in input directory."""
@@ -162,35 +194,76 @@ class LidarArchaeoPipeline:
                 return False
         return True
 
-    def generate_all_visualizations(self, dtm_file, basename, resolution=None):
+    @staticmethod
+    def _expected_output_path(name, basename, file_vis_dir, output_format='avif'):
+        """Return the expected output filename for a visualization step."""
+        ext = 'avif' if output_format == 'avif' else 'webp'
+        if name == 'pos_open':
+            return file_vis_dir / f"{basename}_positive_openness.{ext}"
+        elif name == 'neg_open':
+            return file_vis_dir / f"{basename}_negative_openness.{ext}"
+        elif name == 'hillshade':
+            return file_vis_dir / f"{basename}_hillshade_multi.{ext}"
+        else:
+            return file_vis_dir / f"{basename}_{name}.{ext}"
+
+    def generate_all_visualizations(self, dtm_file, basename, resolution=None, vis_dir=None):
         """Generate all archaeological visualizations for one DTM file.
 
-        Args:
-            resolution: Actual resolution from DTM geotransform. If None, uses self.resolution.
+        Optimisation: SharedDEM is only computed if at least one visualization
+        needs to be generated. When all WebP outputs exist, SharedDEM is
+        skipped entirely (saves ~2min per file on re-runs).
         """
         if resolution is None:
             resolution = self.resolution
         logger.info("  Génération visualisations:")
 
-        # Create per-file subdirectory
-        file_vis_dir = self.vis_dir / basename
+        # Use provided vis_dir (for multi-resolution subdirectories) or default
+        file_vis_dir = vis_dir if vis_dir else (self.vis_dir / basename)
         file_vis_dir.mkdir(exist_ok=True)
+        total = len(self.viz_steps)
 
-        # Pre-compute shared DEM data (gradient, NaN mask, LRM) once for all visualizations
-        logger.info("    Pré-calcul données partagées (gradient, LRM)...")
-        t_shared = time.time()
-        shared = SharedDEM(dtm_file, resolution)
-        logger.info(f"    ✓ Données partagées prêtes ({time.time()-t_shared:.1f}s)")
+        # Phase 1: determine which visualizations need generation
+        needs_generation = {}  # name -> True/False
+        for name, func in self.viz_steps:
+            if self.force:
+                needs_generation[name] = True
+            else:
+                expected_webp = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
+                needs_generation[name] = not expected_webp.exists()
 
+        to_generate = [n for n, needed in needs_generation.items() if needed]
+        ign_only = all(name in ('ortho', 'topo') for name in to_generate)
+        needs_shared = any(name not in ('ortho', 'topo') for name in to_generate)
+
+        if not to_generate:
+            logger.info("    Toutes les visualisations déjà existantes — ignorées")
+            # Still need to return results dict for PDF check
+            vis_results = {}
+            for name, func in self.viz_steps:
+                vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
+            return vis_results
+
+        # Phase 2: compute SharedDEM only if needed
+        shared = None
+        if needs_shared:
+            logger.info("    Pré-calcul données partagées (gradient, LRM)...")
+            t_shared = time.time()
+            shared = SharedDEM(dtm_file, resolution)
+            logger.info(f"    ✓ Données partagées prêtes ({time.time()-t_shared:.1f}s)")
+
+        # Phase 3: generate visualizations
         vis_results = {}
-        total = len(VIZ_STEPS)
+        for idx, (name, func) in enumerate(self.viz_steps, 1):
+            if not needs_generation[name]:
+                logger.info(f"    [{idx}/{total}] {name}: déjà existant, ignoré")
+                vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
+                continue
 
-        for idx, (name, func) in enumerate(VIZ_STEPS, 1):
             # When --force, delete existing TIF to ensure clean regeneration
             if self.force:
                 for tif in file_vis_dir.glob(f"{basename}_{name}.tif"):
                     tif.unlink(missing_ok=True)
-                # Special cases for differently-named TIFs
                 if name == 'pos_open':
                     for tif in file_vis_dir.glob(f"{basename}_positive_openness.tif"):
                         tif.unlink(missing_ok=True)
@@ -201,26 +274,6 @@ class LidarArchaeoPipeline:
                     for tif in file_vis_dir.glob(f"{basename}_hillshade_multi.tif"):
                         tif.unlink(missing_ok=True)
 
-            # Check if output WebP already exists (skip unless --force)
-            if not self.force:
-                # Determine expected WebP filename from the viz name
-                # Special cases for openness and IGN overlays
-                if name == 'pos_open':
-                    expected_webp = file_vis_dir / f"{basename}_positive_openness.webp"
-                elif name == 'neg_open':
-                    expected_webp = file_vis_dir / f"{basename}_negative_openness.webp"
-                elif name == 'hillshade':
-                    expected_webp = file_vis_dir / f"{basename}_hillshade_multi.webp"
-                elif name in ('ortho', 'topo'):
-                    expected_webp = file_vis_dir / f"{basename}_{name}.webp"
-                else:
-                    expected_webp = file_vis_dir / f"{basename}_{name}.webp"
-
-                if expected_webp.exists():
-                    logger.info(f"    [{idx}/{total}] {name}: déjà existant, ignoré")
-                    vis_results[name] = expected_webp  # Track as existing file
-                    continue
-
             logger.info(f"    [{idx}/{total}] {name}...")
             t0 = time.time()
             try:
@@ -242,8 +295,9 @@ class LidarArchaeoPipeline:
             # Free GPU memory between visualizations to prevent OOM
             gpu_cleanup()
 
-        # Convert to WebP (only newly generated TIFs, not skipped ones)
-        logger.info("  Conversion images WebP:")
+        # Convert to output format (only newly generated TIFs, not skipped ones)
+        fmt_label = self.output_format.upper()
+        logger.info(f"  Conversion images {fmt_label}:")
         source_info = {
             'method': self.ground_method,
             'date': datetime.now().strftime('%Y-%m-%d'),
@@ -251,9 +305,9 @@ class LidarArchaeoPipeline:
         }
         for name, tif_file in vis_results.items():
             if tif_file and isinstance(tif_file, Path) and tif_file.suffix == '.tif' and tif_file.exists():
-                webp_file = tif_to_png(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, source_info=source_info)
-                if webp_file:
-                    logger.info(f"    ✓ {webp_file.name}")
+                img_file = tif_to_png(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, source_info=source_info, quality=self.quality, output_format=self.output_format)
+                if img_file:
+                    logger.info(f"    ✓ {img_file.name}")
 
         # Clean up remaining TIF files unless --keep-tif
         if not self.keep_tif:
@@ -262,8 +316,22 @@ class LidarArchaeoPipeline:
 
         return vis_results
 
+    @staticmethod
+    def _res_suffix(resolution):
+        """Return suffix for additional resolutions (empty string for primary)."""
+        if resolution == 0.5:
+            return ""  # Default resolution — no suffix
+        res_str = f"{resolution}".replace('.', 'p')
+        return f"_r{res_str}"
+
     def process_file(self, laz_file):
-        """Process a single LAZ file through the full pipeline."""
+        """Process a single LAZ file through the full pipeline.
+
+        If self.resolutions has multiple entries, processes each resolution:
+        - Primary resolution uses current naming (no suffix)
+        - Additional resolutions use _r0p2 suffix in directories/filenames
+        - Ground classification is done once and shared across resolutions
+        """
         basename = _file_basename(laz_file)
         _file_filter.basename = basename
         t_start = time.time()
@@ -277,74 +345,84 @@ class LidarArchaeoPipeline:
         if not validate_laz(laz_file):
             return False
 
-        # Skip ground classification + DTM if DTM already exists with matching resolution
-        # --force only affects visualizations/PDF, not classification/DTM
-        # Use --force-classification to force reclassification
-        dtm_path = self.dtm_dir / f"{basename}_dtm.tif"
-        if dtm_path.exists():
-            # Check that existing DTM resolution matches requested resolution
-            import rasterio
-            try:
-                with rasterio.open(dtm_path) as src:
-                    existing_res = abs(src.transform.a)
-                if abs(existing_res - self.resolution) > 0.01:
-                    logger.info(f"[1/5] DTM existant à {existing_res}m/px — résolution demandée {self.resolution}m/px → régénération")
+        # Step 1: Ground classification (shared across all resolutions)
+        las_file = None
+        t_classif = 0
+        for i, res in enumerate(self.resolutions):
+            res_suffix = self._res_suffix(res)
+            dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif"
+            if dtm_path.exists():
+                import rasterio
+                try:
+                    with rasterio.open(dtm_path) as src:
+                        existing_res = abs(src.transform.a)
+                    if abs(existing_res - res) > 0.01:
+                        logger.info(f"  DTM{res_suffix} existant à {existing_res}m/px — résolution demandée {res}m/px → régénération")
+                        dtm_path.unlink()
+                    else:
+                        if i == 0:
+                            logger.info(f"[1/5] Classification du sol — sautée (DTM existant)")
+                            logger.info(f"[2/5] Génération DTM {res}m/px — sautée (DTM existant)")
+                        else:
+                            logger.info(f"  DTM {res}m/px déjà existant — ignoré")
+                        continue
+                except Exception:
+                    logger.warning(f"Impossible de lire le DTM existant — régénération")
                     dtm_path.unlink()
-                else:
-                    logger.info(f"[1/5] Classification du sol — sautée (DTM existant à {existing_res}m/px)")
-                    logger.info("[2/5] Génération DTM — sautée (DTM existant)")
-                    dtm_file = dtm_path
-                    t_classif = 0
-                    t_dtm = 0
-            except Exception:
-                logger.warning(f"Impossible de lire le DTM existant — régénération")
-                dtm_path.unlink()
 
-        if not dtm_path.exists():
-            # Step 1: Ground classification
-            logger.info("[1/5] Classification du sol...")
-            t1 = time.time()
-            las_file = classify_ground(laz_file, self.temp_dir, method=self.ground_method, force=self.force_classify)
-            t_classif = time.time() - t1
-            if not las_file:
-                logger.error(f"  ✗ Échec classification ({t_classif:.1f}s)")
-                return False
-            logger.info(f"  ✓ Classification terminée ({t_classif:.1f}s)")
+            # Need to classify/generate DTM for this resolution
+            if las_file is None:
+                # First time: do ground classification
+                logger.info("[1/5] Classification du sol...")
+                t1 = time.time()
+                las_file = classify_ground(laz_file, self.temp_dir, method=self.ground_method, force=self.force_classify)
+                t_classif = time.time() - t1
+                if not las_file:
+                    logger.error(f"  ✗ Échec classification ({t_classif:.1f}s)")
+                    return False
+                logger.info(f"  ✓ Classification terminée ({t_classif:.1f}s)")
 
-            # Step 2: Generate DTM
-            logger.info("[2/5] Génération DTM...")
+            # Generate DTM at this resolution
+            logger.info(f"{'[2/5]' if i == 0 else '  '} Génération DTM {res}m/px...")
             t2 = time.time()
-            dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, self.resolution)
+            dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, res, force=self.force, output_suffix=res_suffix)
             t_dtm = time.time() - t2
             if not dtm_file:
-                logger.error(f"  ✗ Échec DTM ({t_dtm:.1f}s)")
-                return False
-            logger.info(f"  ✓ DTM terminé ({t_dtm:.1f}s)")
+                logger.error(f"  ✗ Échec DTM {res}m/px ({t_dtm:.1f}s)")
+                if i == 0:
+                    return False  # Primary resolution failure is fatal
+                continue  # Additional resolution failure is non-fatal
+            logger.info(f"  ✓ DTM {res}m/px terminé ({t_dtm:.1f}s)")
 
-        # Step 3: Visualizations — use actual resolution from DTM
-        import rasterio
-        with rasterio.open(dtm_file) as src:
-            actual_res = abs(src.transform.a)
-        if abs(actual_res - self.resolution) > 0.01:
-            logger.info(f"  Résolution DTM: {actual_res}m/px (demandée: {self.resolution}m/px)")
-        self.generate_all_visualizations(dtm_file, basename, actual_res)
+        # Process each resolution: visualizations + PDF
+        all_vis_results = {}
+        for res in self.resolutions:
+            res_suffix = self._res_suffix(res)
+            dtm_path = self.dtm_dir / f"{basename}_dtm{res_suffix}.tif"
 
-        # Step 4: PDF report
-        t_pdf = 0
-        file_vis_dir = self.vis_dir / basename
-        pdf_file = self.pdf_dir / f"{basename}_rapport.pdf"
-        if pdf_file.exists() and not self.force:
-            logger.info(f"[4/5] Rapport PDF déjà existant — ignoré: {pdf_file.name}")
-        else:
-            logger.info("[4/5] Rapport PDF A3...")
-            t4 = time.time()
-            generate_pdf_report(basename, file_vis_dir, self.pdf_dir, actual_res)
-            t_pdf = time.time() - t4
-            logger.info(f"  ✓ Rapport PDF terminé ({t_pdf:.1f}s)")
+            if not dtm_path.exists():
+                logger.warning(f"  DTM {res}m/px manquant — visualisations ignorées")
+                continue
+
+            import rasterio
+            with rasterio.open(dtm_path) as src:
+                actual_res = abs(src.transform.a)
+
+            if len(self.resolutions) > 1:
+                logger.info(f"  --- Résolution {res}m/px ---")
+
+            # For additional resolutions, use suffixed subdirectory
+            if res_suffix:
+                vis_dir = self.vis_dir / f"{basename}{res_suffix}"
+            else:
+                vis_dir = self.vis_dir / basename
+
+            vis_dir.mkdir(exist_ok=True)
+
+            self.generate_all_visualizations(dtm_path, basename, actual_res, vis_dir=vis_dir)
 
         t_total = time.time() - t_start
         logger.info(f"✓ {basename} terminé en {t_total:.1f}s")
-        logger.debug(f"  Détails: classification={t_classif:.1f}s, DTM={t_dtm:.1f}s, PDF={t_pdf:.1f}s")
         _file_filter.basename = None
         return True
 
@@ -372,29 +450,42 @@ class LidarArchaeoPipeline:
             logger.info(f"Traitement parallèle avec {self.workers} workers...")
             logger.info(f"Fichiers: {len(files)}")
             with ProcessPoolExecutor(max_workers=self.workers) as executor:
+                # Pass resolutions as comma-separated string for multiprocessing serialization
+                resolutions_str = ','.join(str(r) for r in self.resolutions)
                 future_to_file = {
-                    executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution, self.force, self.ground_method, self.force_classify, self.keep_tif): laz_file
+                    executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), resolutions_str, self.force, self.ground_method, self.force_classify, self.keep_tif, self.quality, self.only_viz, self.skip_viz, self.output_format): laz_file
                     for laz_file in files
                 }
                 done = 0
-                for future in as_completed(future_to_file):
-                    laz_file = future_to_file[future]
-                    done += 1
-                    try:
-                        success = future.result()
-                        results[laz_file.name] = success
-                        status = "✓" if success else "✗"
-                        logger.info(f"  [{done}/{len(files)}] {status} {laz_file.name}")
-                    except Exception as e:
-                        logger.error(f"  [{done}/{len(files)}] ✗ {laz_file.name}: {e}")
-                        logger.debug(f"  Traceback:", exc_info=True)
-                        results[laz_file.name] = False
+                try:
+                    for future in as_completed(future_to_file):
+                        laz_file = future_to_file[future]
+                        done += 1
+                        try:
+                            success = future.result()
+                            results[laz_file.name] = success
+                            status = "✓" if success else "✗"
+                            logger.info(f"  [{done}/{len(files)}] {status} {laz_file.name}")
+                        except Exception as e:
+                            logger.error(f"  [{done}/{len(files)}] ✗ {laz_file.name}: {e}")
+                            logger.debug(f"  Traceback:", exc_info=True)
+                            results[laz_file.name] = False
+                except KeyboardInterrupt:
+                    logger.info("Interruption — annulation des travaux en cours...")
+                    for f in future_to_file:
+                        f.cancel()
+                    executor.shutdown(wait=False, cancel_futures=True)
+                    logger.info("Travaux annulés.")
+                    return
         else:
             total = len(files)
             for idx, laz_file in enumerate(files, 1):
                 logger.info(f"--- Fichier {idx}/{total} ---")
                 try:
                     results[laz_file.name] = self.process_file(laz_file)
+                except KeyboardInterrupt:
+                    logger.info("Interruption — arrêt immédiat.")
+                    return
                 except Exception as e:
                     logger.error(f"✗ Erreur traitement {laz_file.name}: {e}")
                     logger.debug("Traceback:", exc_info=True)
@@ -408,18 +499,18 @@ class LidarArchaeoPipeline:
         logger.info("=" * 60)
         logger.info("RÉSUMÉ")
         logger.info("=" * 60)
+        for name, ok in results.items():
+            status = "✓" if ok else "✗"
+            logger.info(f"  {status} {name}")
+        logger.info("-" * 60)
         logger.info(f"  Succès : {success_count}/{len(results)}")
         if fail_count:
             logger.info(f"  Échecs : {fail_count}/{len(results)}")
-            for name, ok in results.items():
-                if not ok:
-                    logger.info(f"    ✗ {name}")
         logger.info(f"  Durée totale : {t_pipeline_total:.1f}s ({t_pipeline_total/60:.1f}min)")
 
         logger.info(f"\nRésultats dans: {self.output_dir}")
         logger.info(f"  • DTM           : {self.dtm_dir}")
         logger.info(f"  • Visualisations: {self.vis_dir}")
-        logger.info(f"  • Rapports PDF  : {self.pdf_dir}")
 
         # Clean up temporary files
         logger.info("Nettoyage des fichiers temporaires...")
@@ -435,7 +526,7 @@ class LidarArchaeoPipeline:
             logger.warning(f"  Note: Impossible de supprimer les fichiers temporaires: {e}")
 
 
-def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False):
+def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False, quality=98, only_viz=None, skip_viz=None, output_format='avif'):
     """Standalone function for multiprocessing — creates its own pipeline instance.
 
     Each worker gets its own temp directory to avoid file conflicts.
@@ -456,7 +547,7 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo
         worker_logger.addHandler(handler)
     worker_logger.addFilter(_file_filter)
 
-    pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif)
+    pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format)
     basename = _file_basename(laz_file_str)
     pipeline.temp_dir = pipeline.output_dir / "temp" / basename
     pipeline.temp_dir.mkdir(exist_ok=True)

lidar_pipeline/rendering.py

@@ -1,8 +1,8 @@
-"""Rendering module: colormap registry, GeoTIFF-to-WebP conversion, and PDF report generation.
+"""Rendering module: colormap registry, GeoTIFF-to-image conversion, and PDF report generation.
 
 Contains:
 - COLORMAPS: registry mapping filename keywords to (cmap, title, legend, description)
-- tif_to_png(): convert a GeoTIFF to a WebP visualization with legend, scale bar, north arrow
+- tif_to_png(): convert a GeoTIFF to a WebP/AVIF visualization with legend, scale bar, north arrow
 - generate_pdf_report(): generate an A3 PDF report with all visualizations
 """
 
@@ -74,18 +74,10 @@ COLORMAPS = {
         'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées',
         'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
     },
-    'svf': {
-        'cmap': 'viridis',
-        'title': 'Sky-View Factor (Ray-tracing 16 directions)',
-        'legend': 'Fraction de ciel visible depuis chaque point\nSombre = Encaissé (fossés, vallées, rues)\nClair = Dégagé (sommets, plateformes, plateaux)',
-        'description': 'Ray-tracing sur 16 azimuts — élimine l\'éclairage, détecte structures linéaires et enclos',
-        'vmin_mode': 'percentile', 'vmin_pct': 5,
-        'vmax_mode': 'percentile', 'vmax_pct': 95,
-    },
     'mslrm': {
         'cmap': 'RdBu_r',
-        'title': 'MSRM - Multi-Scale Relief Model (5 échelles)',
-        'legend': 'Relief combiné à 5 échelles (5m à 100m)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve, fossé)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = 5 échelles combinées\nMSRM détecte du micro au macro',
+        'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)',
+        'legend': 'Relief combiné multi-échelles (adapté à la résolution)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nMSRM détecte du micro au macro',
         'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
@@ -114,8 +106,8 @@ COLORMAPS = {
     },
     'tpi': {
         'cmap': 'BrBG',
-        'title': 'TPI - Topographic Position Index (2 échelles)',
-        'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\nCombine échelle fine 5m + large 100m',
+        'title': 'TPI - Topographic Position Index (4 échelles)',
+        'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\nCombine 4 échelles : 3m, 15m, 50m, 200m',
         'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
@@ -128,23 +120,23 @@ COLORMAPS = {
     },
     'roughness': {
         'cmap': 'magma',
-        'title': 'Rugosité de Surface (Écart-type local 5m)',
-        'legend': 'Irrégularité du terrain dans un voisinage de 5m\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nMax: {vmax:.2f}m',
+        'title': 'Rugosité Multi-Échelle (3m + 15m)',
+        'legend': 'Irrégularité du terrain combinée fine + large\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nCombine rugosité fine 3m (70%) + large 15m (30%)',
         'description': 'Mesure la variabilité locale — surfaces anthropiques lisses vs naturelles rugueuses',
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 97,
     },
     'anomalies': {
         'cmap': 'coolwarm',
-        'title': 'Anomalies Statistiques (Z-score x Moran\'s I)',
-        'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine z-score (intensité) et\nMoran\'s I (regroupement spatial)',
+        'title': 'Anomalies Statistiques (MSRM multi-échelle + Moran\'s I)',
+        'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine MSRM normalisé (intensité) et\nMoran\'s I (regroupement spatial)',
         'description': 'Détecte uniquement les anomalies statistiquement significatives — filtre le bruit de fond',
         'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
     },
     'wavelet': {
         'cmap': 'cividis',
         'title': 'Ondelette Mexican Hat (CWT multi-échelle)',
-        'legend': 'Réponse de la transformée en ondelette à 5 échelles\nÉchelles: 2m, 5m, 10m, 20m, 50m\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
+        'legend': 'Réponse de la transformée en ondelette\nÉchelles adaptées à la résolution\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
         'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
@@ -156,14 +148,6 @@ COLORMAPS = {
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 98,
     },
-    'local_dominance': {
-        'cmap': 'RdYlBu_r',
-        'title': 'Dominance Locale (position relative dans le voisinage)',
-        'legend': 'Proportion du voisinage sous le point central\nRouge = Point dominant (sommet, crête)\nBleu = Point encaissé (fossé, vallée)\nRayon: 15m',
-        'description': 'Mesure la saillie locale — complémentaire de l\'openness',
-        'vmin_mode': 'percentile', 'vmin_pct': 2,
-        'vmax_mode': 'percentile', 'vmax_pct': 98,
-    },
 }
 
 # RGB entries (ortho/topo) are handled specially
@@ -271,22 +255,26 @@ def _nice_scale(extent_m):
     return chosen, f"{chosen} m"
 
 
-def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None):
-    """Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar.
+def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, quality=98, output_format='avif'):
+    """Convert GeoTIFF to visualization image (WebP or AVIF) with GPS coordinates, legend, and scale bar.
 
     Args:
         tif_file: Path to input GeoTIFF.
-        vis_dir: Output directory for the WebP file.
+        vis_dir: Output directory for the image file.
         resolution: Grid resolution in m/px.
         keep_tif: If True, keep the source TIFF after conversion.
+        source_info: Dict with method/date/basename for metadata.
+        quality: Image quality (1-100). Use 100 for lossless. Default 95.
+        output_format: Output format ('webp' or 'avif'). Default 'webp'.
 
     Returns:
-        Path to output WebP file, or None on failure.
+        Path to output image file, or None on failure.
     """
     if not tif_file or not tif_file.exists():
         return None
 
-    webp_file = vis_dir / f"{tif_file.stem}.webp"
+    ext = 'avif' if output_format == 'avif' else 'webp'
+    output_file = vis_dir / f"{tif_file.stem}.{ext}"
 
     try:
         with rasterio.open(tif_file) as src:
@@ -580,7 +568,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None):
 
             fig.patch.set_facecolor('white')
 
-            # Save as PNG then convert to WebP — fixed layout, no bbox_inches='tight'
+            # Save as PNG then convert to final format — fixed layout, no bbox_inches='tight'
             save_dpi = 200 if width > 3000 else 150
             png_temp = vis_dir / f"{tif_file.stem}_temp.png"
             try:
@@ -589,17 +577,21 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None):
                 plt.close()
 
             img = PILImage.open(str(png_temp))
-            img.save(str(webp_file), format='WEBP', lossless=True)
+            pil_format = 'AVIF' if output_format == 'avif' else 'WEBP'
+            if quality >= 100:
+                img.save(str(output_file), format=pil_format, lossless=True)
+            else:
+                img.save(str(output_file), format=pil_format, quality=quality)
             png_temp.unlink(missing_ok=True)
 
             # Delete source TIFF (unless --keep-tif)
             if not keep_tif:
                 tif_file.unlink(missing_ok=True)
 
-            return webp_file
+            return output_file
 
     except Exception as e:
-        logger.error(f"    Erreur conversion WebP: {e}", exc_info=True)
+        logger.error(f"    Erreur conversion {ext.upper()}: {e}", exc_info=True)
         return None
 
 

lidar_pipeline/visualizations.py

@@ -33,10 +33,13 @@ else:
 class SharedDEM:
     """Pre-computed DEM data shared across all visualizations.
 
-    Reads the DEM once and pre-computes:
+    Reads the DEM once and lazily computes on first access:
     - NaN mask and filled DEM (avoids 20+ calls to _fill_nans)
     - Gradient components (shared by hillshade, slope, aspect, curvature)
     - LRM at 15m kernel (shared by lrm + anomalies)
+
+    Attributes are computed lazily on first access to avoid computing
+    data that is never used (e.g. LRM when only hillshade needs generation).
     """
 
     def __init__(self, dem_file, resolution):
@@ -48,25 +51,69 @@ class SharedDEM:
         self.nan_mask = np.isnan(dem_np)
         self.dem_np = dem_np.astype(np.float32)
 
-        # Pre-fill NaNs once (saves ~20 calls to NearestNDInterpolator)
-        self.filled, _ = _fill_nans(self.dem_np)
+        # Lazy caches — computed on first access
+        self._filled = None
+        self._gradient = None  # (dy, dx, slope_rad, slope_deg, aspect)
+        self._lrm_15 = None
 
-        # Initialize GPU lazy caches before any filter calls
+        # GPU lazy caches
         self._filled_gpu = None
         self._dem_gpu = None
 
-        # Pre-compute gradient (shared by hillshade, slope, aspect, curvature)
-        self.dy = np.gradient(self.filled, resolution, axis=0)
-        self.dx = np.gradient(self.filled, resolution, axis=1)
-        self.slope_rad = np.arctan(np.sqrt(self.dx**2 + self.dy**2))
-        self.slope_deg = np.degrees(self.slope_rad)
-        self.aspect = np.mod(np.degrees(np.arctan2(self.dy, self.dx)), 360)
+    @property
+    def filled(self):
+        """Filled DEM (NaN interpolated) — computed lazily."""
+        if self._filled is None:
+            logger.debug("    → Calcul filled DEM (interpolation NaN)...")
+            self._filled, _ = _fill_nans(self.dem_np)
+        return self._filled
 
-        # Pre-compute LRM at 15m (shared by lrm + anomalies)
-        sigma_15 = 15.0 / resolution
-        local_mean_15 = _filter_nanaware_from_filled(self, xp_gaussian_filter, sigma=sigma_15)
-        self.lrm_15 = self.dem_np - local_mean_15
-        self.lrm_15[self.nan_mask] = np.nan
+    @property
+    def dy(self):
+        self._ensure_gradient()
+        return self._gradient[0]
+
+    @property
+    def dx(self):
+        self._ensure_gradient()
+        return self._gradient[1]
+
+    @property
+    def slope_rad(self):
+        self._ensure_gradient()
+        return self._gradient[2]
+
+    @property
+    def slope_deg(self):
+        self._ensure_gradient()
+        return self._gradient[3]
+
+    @property
+    def aspect(self):
+        self._ensure_gradient()
+        return self._gradient[4]
+
+    @property
+    def lrm_15(self):
+        """LRM at 15m kernel — computed lazily."""
+        if self._lrm_15 is None:
+            logger.debug("    → Calcul LRM 15m...")
+            sigma_15 = 15.0 / self.resolution
+            local_mean_15 = _filter_nanaware_from_filled(self, xp_gaussian_filter, sigma=sigma_15)
+            self._lrm_15 = self.dem_np - local_mean_15
+            self._lrm_15[self.nan_mask] = np.nan
+        return self._lrm_15
+
+    def _ensure_gradient(self):
+        """Compute gradient components lazily on first access."""
+        if self._gradient is None:
+            logger.debug("    → Calcul gradient...")
+            dy = np.gradient(self.filled, self.resolution, axis=0)
+            dx = np.gradient(self.filled, self.resolution, axis=1)
+            slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
+            slope_deg = np.degrees(slope_rad)
+            aspect = np.mod(np.degrees(np.arctan2(dy, dx)), 360)
+            self._gradient = (dy, dx, slope_rad, slope_deg, aspect)
 
     @property
     def filled_gpu(self):
@@ -203,7 +250,7 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
 def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
     """Generate multi-directional hillshade with contrast enhancement — GPU if available.
 
-    Combines 4-direction hillshade (NW, NE, SW, SE) with slope shading.
+    Combines 8-direction hillshade with slope shading for balanced illumination.
     Applies percentile normalization and gamma correction to restore
     contrast lost by averaging multiple azimuths.
     """
@@ -232,8 +279,9 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
             sin_slope = xp.sin(slope)
             cos_slope = xp.cos(slope)
 
-        azimuts = [315, 45, 225, 135]
-        altitude = 30
+        # 8 azimuths for balanced illumination (eliminates directional bias)
+        azimuts = [0, 45, 90, 135, 180, 225, 270, 315]
+        altitude = 35  # Higher altitude for better micro-relief detection
         hillshades = []
 
         alt_rad = xp.radians(xp.array(altitude))
@@ -250,7 +298,6 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
         combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
 
         # Contrast enhancement: percentile stretch + gamma
-        # Averaging 4 azimuths flattens contrast — this restores it
         combined_np = to_cpu(combined)
         nan_mask = shared.nan_mask if shared else np.isnan(to_cpu(dem_np) if HAS_GPU else dem_np)
         valid = combined_np[~nan_mask]
@@ -368,7 +415,11 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None):
 # ============================================================
 
 def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None):
-    """Local Relief Model - deviation from local mean (GPU if available)."""
+    """Local Relief Model - deviation from local mean (GPU if available).
+
+    Kernel sigma adapts to resolution: finer kernel at higher resolution
+    to capture micro-relief details. At 0.5m/px: 15m, at 0.2m/px: ~5m.
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Local Relief Model{gpu_tag}...")
     t0 = time.time()
@@ -382,7 +433,10 @@ def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None):
         else:
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
-            local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15/resolution)
+            # Adapt sigma to resolution: standard 15m at 0.5m, finer at higher res
+            sigma_m = max(5.0, 15.0 * 0.5 / resolution)
+            logger.info(f"    LRM sigma={sigma_m:.1f}m (résolution {resolution}m/px)")
+            local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_m / resolution)
             lrm = dem_np - local_mean
             lrm[nan_mask] = np.nan
         _save_tif(output, lrm.astype(np.float32), transform, crs)
@@ -426,7 +480,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
         angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
         dx_dir = np.cos(angles)
         dy_dir = np.sin(angles)
-        max_dist = int(50 / res)
+        max_dist = int(100 / res)
 
         padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
         svf = xp.zeros_like(dem)
@@ -473,6 +527,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
     - 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.
+    Ray radius adapts to resolution: 100m for better detection of large enclosures.
     """
     name = "positive_openness" if positive else "negative_openness"
     gpu_tag = " [GPU]" if HAS_GPU else ""
@@ -501,7 +556,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
         angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
         dx_dir = np.cos(angles)
         dy_dir = np.sin(angles)
-        max_dist = int(50 / res)
+        max_dist = int(100 / res)
 
         padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
         openness_sum = xp.zeros_like(dem)
@@ -599,7 +654,11 @@ def generate_local_dominance(dem_file, basename, vis_dir, resolution, shared=Non
 
 
 def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
-    """Multi-Scale Relief Model (MSRM) - LRM at 5 scales combined (GPU if available)."""
+    """Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).
+
+    Scales adapt to resolution. Std normalization per scale.
+    Weighted combination favoring archaeologically relevant scales (5-25m).
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Multi-Scale Relief Model (MSRM){gpu_tag}...")
     t0 = time.time()
@@ -615,7 +674,18 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
 
-        sigmas = [5, 10, 25, 50, 100]
+        # Adaptive scales: finer at higher resolution
+        min_scale = max(2.0, resolution * 4)
+        candidate_scales = [2, 5, 10, 20, 50, 100, 200]
+        sigmas = [s for s in candidate_scales if s >= min_scale]
+
+        # Archaeological weights: favor 5-25m range (ditches, enclosures, tumulus)
+        scale_weights = {
+            2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6, 200: 0.4,
+        }
+        weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
+
+        logger.info(f"    MSRM échelles: {sigmas}m")
         lrm_stack = []
 
         for sigma in sigmas:
@@ -626,16 +696,19 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
                 local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
             lrm = dem_np - local_mean
             lrm[nan_mask] = np.nan
+            # Std normalization: x / std — preserves sign and contrast better than z-score
             valid_lrm = lrm[~nan_mask]
             lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
-            lrm_norm = lrm / lrm_std
-            lrm_stack.append(lrm_norm.astype(np.float32))
+            lrm = lrm / lrm_std
+            lrm_stack.append(lrm.astype(np.float32))
 
+        # Weighted combination
         lrm_array = np.array(lrm_stack)
+        weights_3d = weights[:, np.newaxis, np.newaxis]
         with np.errstate(invalid='ignore', divide='ignore'):
             with warnings.catch_warnings():
                 warnings.filterwarnings('ignore', message='Mean of empty slice')
-                mslrm = np.sqrt(np.nanmean(lrm_array ** 2, axis=0))
+                mslrm = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
         mslrm[nan_mask] = np.nan
         _save_tif(output, mslrm.astype(np.float32), transform, crs)
         logger.info(f"  ✓ MSRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
@@ -649,7 +722,8 @@ def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None):
     """Multi-Scale Topographic Position Index (GPU if available).
 
     TPI = elevation - mean(neighborhood).
-    Computed at fine (5m) and broad (100m) scales.
+    Computed at 4 scales with std normalization and weighted combination.
+    Weights favor fine and medium scales (archaeologically relevant).
     """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → TPI multi-échelle{gpu_tag}...")
@@ -666,29 +740,34 @@ def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None):
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
 
-        fine_size = int(5 / resolution)
-        if fine_size % 2 == 0:
-            fine_size += 1
-        if shared:
-            fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
-        else:
-            fine_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=fine_size)
-        tpi_fine = dem_np - fine_mean
-        tpi_fine[nan_mask] = np.nan
+        # 4 scales: fine (3m), medium (15m), broad (50m), landscape (200m)
+        scales_m = [3, 15, 50, 200]
+        weights = [1.5, 2.0, 1.2, 0.5]  # Favor medium scales (ditches, enclosures)
 
-        broad_size = int(100 / resolution)
-        if broad_size % 2 == 0:
-            broad_size += 1
-        if shared:
-            broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size)
-        else:
-            broad_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=broad_size)
-        tpi_broad = dem_np - broad_mean
-        tpi_broad[nan_mask] = np.nan
+        tpi_stack = []
+        for scale_m, weight in zip(scales_m, weights):
+            size = max(3, int(scale_m / resolution))
+            if size % 2 == 0:
+                size += 1
+            if shared:
+                local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=size)
+            else:
+                local_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=size)
+            tpi = dem_np - local_mean
+            tpi[nan_mask] = np.nan
+            # Std normalization — preserves sign and contrast better than z-score
+            valid = tpi[~nan_mask]
+            tpi_std = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
+            tpi = tpi / tpi_std
+            tpi_stack.append(tpi.astype(np.float32))
 
-        fine_std = max(np.nanstd(tpi_fine[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
-        broad_std = max(np.nanstd(tpi_broad[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
-        tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std)
+        # Weighted combination
+        tpi_array = np.array(tpi_stack)
+        weights_3d = np.array(weights)[:, np.newaxis, np.newaxis]
+        with np.errstate(invalid='ignore', divide='ignore'):
+            with warnings.catch_warnings():
+                warnings.filterwarnings('ignore', message='Mean of empty slice')
+                tpi_combined = np.nansum(tpi_array * weights_3d, axis=0) / np.sum(weights)
         tpi_combined[nan_mask] = np.nan
 
         _save_tif(output, tpi_combined.astype(np.float32), transform, crs)
@@ -709,7 +788,7 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
     """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.
+    steep areas get smaller kernels. Scales adapt to resolution.
     """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → SAILORE (LRM adaptatif){gpu_tag}...")
@@ -731,8 +810,13 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
             slope_deg = np.degrees(slope)
             slope_deg[nan_mask] = np.nan
 
-        sigma_min = 2.0 / resolution
-        sigma_max = 25.0 / resolution
+        # Adaptive scales: finer at higher resolution
+        sigma_min_m = max(1.0, 2.0 * 0.5 / resolution)  # 2m at 0.5, ~5m at 0.2
+        sigma_mid_m = max(5.0, 13.5 * 0.5 / resolution)  # 13.5m at 0.5, ~33m at 0.2
+        sigma_max_m = max(5.0, 25.0 * 0.5 / resolution)  # 25m at 0.5, ~62m at 0.2
+        sigma_min = sigma_min_m / resolution
+        sigma_max = sigma_max_m / resolution
+        sigma_mid = (sigma_min + sigma_max) / 2
         slope_norm = np.clip(slope_deg / 30.0, 0, 1)
 
         if shared:
@@ -775,7 +859,11 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
 # ============================================================
 
 def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
-    """Surface roughness - standard deviation of elevation in a window (GPU-accelerated)."""
+    """Surface roughness - multi-scale standard deviation (GPU-accelerated).
+
+    Combines fine (3m) and broad (15m) roughness for better detection
+    of archaeological features at multiple scales.
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Rugosité de surface{gpu_tag}...")
     t0 = time.time()
@@ -791,20 +879,43 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
 
-        window_size = int(5 / resolution)
-        if window_size % 2 == 0:
-            window_size += 1
+        # Fine roughness (3m window)
+        fine_size = max(3, int(3 / resolution))
+        if fine_size % 2 == 0:
+            fine_size += 1
 
-        # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (NaN-aware)
         if shared:
-            local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window_size)
-            # For local_mean_sq, we need to filter filled², not filled
-            local_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=window_size)
-            local_mean_sq[shared.nan_mask] = np.nan
+            fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
+            fine_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
+            fine_mean_sq[shared.nan_mask] = np.nan
         else:
-            local_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=window_size)
-            local_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=window_size)
-        roughness = np.sqrt(np.maximum(local_mean_sq - local_mean * local_mean, 0))
+            fine_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=fine_size)
+            fine_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
+        roughness_fine = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0))
+        roughness_fine[nan_mask] = np.nan
+
+        # Broad roughness (15m window)
+        broad_size = max(3, int(15 / resolution))
+        if broad_size % 2 == 0:
+            broad_size += 1
+
+        if shared:
+            broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size)
+            broad_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
+            broad_mean_sq[shared.nan_mask] = np.nan
+        else:
+            broad_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=broad_size)
+            broad_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
+        roughness_broad = np.sqrt(np.maximum(broad_mean_sq - broad_mean * broad_mean, 0))
+        roughness_broad[nan_mask] = np.nan
+
+        # Std normalization per scale then weighted combination
+        fine_valid = roughness_fine[~nan_mask]
+        broad_valid = roughness_broad[~nan_mask]
+        fine_std = max(np.nanstd(fine_valid), 0.01) if len(fine_valid) > 0 else 0.01
+        broad_std = max(np.nanstd(broad_valid), 0.01) if len(broad_valid) > 0 else 0.01
+
+        roughness = 0.7 * roughness_fine / fine_std + 0.3 * roughness_broad / broad_std
         roughness[nan_mask] = np.nan
 
         roughness = to_cpu(roughness)
@@ -821,7 +932,11 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
 # ============================================================
 
 def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
-    """Statistical anomaly detection - z-score of local relief + Local Moran's I — GPU if available."""
+    """Statistical anomaly detection - std-normalized multi-scale relief + Local Moran's I — GPU if available.
+
+    Uses MSRM (multi-scale LRM) instead of single-scale LRM for better detection
+    of anomalies at all scales.
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Détection anomalies statistiques{gpu_tag}...")
     t0 = time.time()
@@ -833,30 +948,60 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
             crs = shared.crs
             dem_np = shared.dem_np
             nan_mask = shared.nan_mask
-            lrm = shared.lrm_15.copy()
         else:
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
-            lrm_mean_val = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15 / resolution)
-            lrm = dem_np - lrm_mean_val
+
+        # Multi-scale LRM: compute MSRM-like combined relief
+        min_scale = max(2.0, resolution * 4)
+        candidate_scales = [2, 5, 10, 20, 50, 100]
+        sigmas = [s for s in candidate_scales if s >= min_scale]
+        lrm_stack = []
+
+        for sigma in sigmas:
+            sigma_px = sigma / resolution
+            if shared:
+                local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px)
+            else:
+                local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
+            lrm = dem_np - local_mean
             lrm[nan_mask] = np.nan
+            # Std normalization — preserves contrast better than z-score
+            valid_lrm = lrm[~nan_mask]
+            lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
+            lrm_norm = lrm / lrm_std
+            else:
+                lrm_norm = lrm
+            lrm_stack.append(lrm_norm.astype(np.float32))
 
-        valid_lrm = lrm[~nan_mask]
-        lrm_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0
-        lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
-        z_score = (lrm - lrm_mean) / lrm_std
+        # Weighted RMS combination (favor 5-25m scales)
+        scale_weights = {2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6}
+        weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
+        lrm_array = np.array(lrm_stack)
+        weights_3d = weights[:, np.newaxis, np.newaxis]
+        with np.errstate(invalid='ignore', divide='ignore'):
+            with warnings.catch_warnings():
+                warnings.filterwarnings('ignore', message='Mean of empty slice')
+                msrm = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
+        msrm[nan_mask] = np.nan
 
-        window = int(10 / resolution)
+        # Std normalization of MSRM — preserves contrast better than z-score
+        valid_msrm = msrm[~nan_mask]
+        msrm_std = max(np.nanstd(valid_msrm), 0.01) if len(valid_msrm) > 0 else 0.01
+        z_score = msrm / msrm_std
+
+        # Local Moran's I for spatial clustering
+        window = max(3, int(10 / resolution))
         if window % 2 == 0:
             window += 1
 
         if shared:
-            local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window)
+            local_mean_z = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window)
         else:
-            local_mean = _filter_nanaware(z_score, xp_uniform_filter, size=window)
-        z_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0
-        z_std = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
-        morans_i = z_score * (local_mean - z_mean) / z_std
+            local_mean_z = _filter_nanaware(z_score, xp_uniform_filter, size=window)
+        z_mean_global = np.nanmean(z_score[~nan_mask]) if np.any(~nan_mask) else 0
+        z_std_global = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
+        morans_i = z_score * (local_mean_z - z_mean_global) / z_std_global
         anomaly_score = np.abs(z_score) * np.sign(morans_i)
         anomaly_score[nan_mask] = np.nan
 
@@ -875,7 +1020,13 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
 def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
     """Mexican Hat wavelet multi-scale analysis (GPU if available).
 
-    CWT 2D at multiple scales to detect circular features.
+    CWT 2D at multiple scales adapted to resolution.
+    - At 0.5m/px: [1, 2, 5, 10, 20, 50, 100]m
+    - At 0.2m/px: [0.5, 1, 2, 5, 10, 20, 50, 100]m
+    - Higher resolution = more fine scales available
+
+    Uses std normalization per scale and weighted combination
+    with emphasis on archaeologically relevant scales (2-50m).
     """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Ondelette Mexican Hat multi-échelle{gpu_tag}...")
@@ -894,7 +1045,25 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
             nan_mask = np.isnan(dem_np)
             filled, _ = _fill_nans(dem_np.astype(np.float64))
 
-        scales = [2, 5, 10, 20, 50]
+        # Adapt scales to resolution: finer scales available at higher resolution
+        min_scale = max(resolution * 2, 1.0)
+        candidate_scales = [0.5, 1, 2, 5, 10, 20, 50, 100]
+        scales = [s for s in candidate_scales if s >= min_scale]
+
+        # Weights favor archaeological scales (2-50m: ditches, enclosures, tumulus)
+        scale_weights = {
+            0.5: 0.6,    # Fine texture
+            1.0: 0.8,    # Micro-relief
+            2.0: 1.5,    # Small ditches, paths — key scale
+            5.0: 2.0,    # Fossés, small enclosures — key archaeological scale
+            10.0: 1.8,   # Medium structures
+            20.0: 1.5,   # Large enclosures, tumulus
+            50.0: 1.0,   # Very large enclosures
+            100.0: 0.6,  # Landscape-level features
+        }
+        weights = np.array([scale_weights.get(s, 1.0) for s in scales])
+
+        logger.info(f"    Échelles CWT: {scales}m (résolution {resolution}m/px)")
         wavelet_stack = []
 
         for scale_m in scales:
@@ -907,15 +1076,21 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
                 from scipy.ndimage import gaussian_laplace
                 response = -gaussian_laplace(filled, sigma=sigma_px)
             response[nan_mask] = np.nan
+
+            # Std normalization: scale by standard deviation to make scales comparable
             valid = response[~nan_mask]
             std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
             response = response / std_val
             wavelet_stack.append(response)
 
+        # Weighted RMS: sqrt(sum(w * x²) / sum(w))
+        # Preserves contrast at key archaeological scales
+        stack = np.array(wavelet_stack)
+        weights_3d = weights[:, np.newaxis, np.newaxis]
         with np.errstate(invalid='ignore', divide='ignore'):
             with warnings.catch_warnings():
                 warnings.filterwarnings('ignore', message='Mean of empty slice')
-                combined = np.sqrt(np.nanmean(np.array(wavelet_stack) ** 2, axis=0))
+                combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(weights))
         combined[nan_mask] = np.nan
 
         _save_tif(output, combined.astype(np.float32), transform, crs)

run.sh

@@ -3,18 +3,21 @@
 # Utilisation: ./run.sh [options]
 #
 # Options:
-#   -r RESOLUTION   Résolution en m/px (défaut: 0.5)
+#   -r RESOLUTION   Résolution en m/px, ou multiples séparées par virgules (défaut: 0.5, ex: 0.5,0.2)
 #   -w WORKERS      Nombre de workers parallèles (défaut: 1)
 #   -g              Activer l'accélération GPU
-#   -v              Mode verbeux (timestamps + niveaux)
-#   --debug         Mode debug (détails internes fichier:ligne)
-#   -f / --force    Régénérer tous les fichiers même si existants
-#   --keep-tif       Conserver les fichiers TIFF pour régénérer les WebP
+#   -v              Mode verbeux
+#   --debug         Mode debug
+#   -f / --force    Régénérer tous les fichiers
+#   --keep-tif      Conserver les fichiers TIFF
 #   --force-classification
-#                   Reclassifier le sol même si le fichier .las existe déjà
 #   --ground-classification {auto,smrf,csf}
-#                   Méthode de classification du sol (défaut: auto)
-#   --file NOM...   Traiter un ou plusieurs fichiers LAZ spécifiques
+#   --quality N     Qualité image 1-100 (défaut: 98)
+#   --lossless      Compression lossless
+#   --format FMT    Format de sortie : avif (défaut) ou webp
+#   --only VIZ...   Générer uniquement ces visualisations
+#   --skip VIZ...   Exclure ces visualisations
+#   --file NOM...   Traiter un ou plusieurs fichiers LAZ
 #   --test          Exécuter les tests unitaires
 #   -h              Afficher l'aide complète
 
@@ -32,7 +35,7 @@ if [ $# -eq 0 ]; then
     echo "Usage: $0 [options]"
     echo ""
     echo "Options:"
-    echo "  -r RESOLUTION   Résolution en m/px (défaut: 0.5)"
+    echo "  -r RESOLUTION   Résolution en m/px, ou multiples (défaut: 0.5, ex: 0.5,0.2)"
     echo "  -w WORKERS      Nombre de workers CPU parallèles (défaut: 1)"
     echo "  -g              Activer l'accélération GPU NVIDIA"
     echo "  -v              Mode verbeux (timestamps + niveaux)"
@@ -52,6 +55,7 @@ if [ $# -eq 0 ]; then
     echo "  $0 -g -w 4                               # GPU + 4 workers"
     echo "  $0 -g -v                                 # GPU + verbeux"
     echo "  $0 -g -r 0.2                             # Haute résolution"
+    echo "  $0 -g -r 0.5,0.2                          # Multi-résolution (0.5m + 0.2m)"
     echo "  $0 -g --force                            # Régénérer WebP (DTM conservé si --keep-tif)"
     echo "  $0 -g --force-classification             # Reclassifier le sol seulement"
     echo "  $0 -g --ground-classification csf        # Forcer CSF (terrain complexe)"
@@ -68,6 +72,10 @@ FILE_ARGS=""
 GROUND_METHOD=""
 FORCE_CLASSIFY_FLAG=""
 KEEP_TIF_FLAG=""
+QUALITY=""
+FORMAT_FLAG=""
+ONLY_FLAG=""
+SKIP_FLAG=""
 
 # Parse arguments manually (more robust than getopts for mixed short/long options)
 while [ $# -gt 0 ]; do
@@ -83,6 +91,11 @@ while [ $# -gt 0 ]; do
         --keep-tif) KEEP_TIF_FLAG="--keep-tif"; shift ;;
         --ground-classification) GROUND_METHOD="$2"; shift 2 ;;
         --ground-classification=*) GROUND_METHOD="${1#--ground-classification=}"; shift ;;
+        --quality) QUALITY="--quality $2"; shift 2 ;;
+        --lossless) QUALITY="--lossless"; shift ;;
+        --format) FORMAT_FLAG="--format $2"; shift 2 ;;
+        --only) shift; ONLY_FLAG="--only"; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do ONLY_FLAG="$ONLY_FLAG $1"; shift; done ;;
+        --skip) shift; SKIP_FLAG="--skip"; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do SKIP_FLAG="$SKIP_FLAG $1"; shift; done ;;
         --file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;;
         --test) ;;  # Handled below
         -h|--help|-help)
@@ -102,17 +115,30 @@ while [ $# -gt 0 ]; do
             echo "  --keep-tif      Conserver les TIFF pour régénérer les WebP"
             echo "  --ground-classification {auto,smrf,csf}"
             echo "                  Méthode de classification du sol (défaut: auto)"
-            echo "  --file NOM...   Traiter un ou plusieurs fichiers LAZ (nom complet sans .laz/.las)"
+            echo "  --quality N     Qualité image 1-100 (défaut: 98, 100=lossless)"
+            echo "  --lossless      Compression lossless (équivalent à --quality 100)"
+            echo "  --format FMT    Format de sortie : avif (défaut) ou webp"
+            echo "  --only VIZ...   Générer uniquement ces visualisations"
+            echo "  --skip VIZ...   Exclure ces visualisations"
+            echo "  --file NOM...   Traiter un ou plusieurs fichiers LAZ"
             echo "  --test          Exécuter les tests unitaires"
             echo "  -h              Afficher cette aide"
             echo ""
+            echo "Visualisations disponibles:"
+            echo "  hillshade slope aspect curvature lrm pos_open neg_open"
+            echo "  mslrm tpi sailore roughness anomalies wavelet flow ortho topo"
+            echo ""
             echo "Exemples:"
             echo "  $0 -g                                    # GPU, auto"
             echo "  $0 -g -w 4                               # GPU + 4 workers"
             echo "  $0 -g -v                                 # GPU + verbeux"
             echo "  $0 -g -r 0.2                             # Haute résolution"
-            echo "  $0 -g --force                            # Régénérer WebP (DTM conservé si --keep-tif)"
-            echo "  $0 -g --force-classification             # Reclassifier le sol seulement"
+    echo "  $0 -g -r 0.5,0.2                          # Multi-résolution (0.5m + 0.2m)"
+            echo "  $0 -g --force                            # Régénérer WebP"
+            echo "  $0 -g --only hillshade svf lrm           # Seulement 3 visualisations"
+            echo "  $0 -g --skip ortho topo                  # Sans les overlays IGN"
+            echo "  $0 -g --lossless                         # WebP lossless"
+            echo "  $0 -g --quality 90                        # WebP qualité 90"
             echo "  $0 -g --ground-classification csf        # Forcer CSF (terrain complexe)"
             echo "  $0 -g --file LHD_...IGN69.copc           # Un fichier"
             exit 0
@@ -160,16 +186,29 @@ echo "  Verbeux    : $([ -n "$VERBOSE_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Force      : $([ -n "$FORCE_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Force classif.: $([ -n "$FORCE_CLASSIFY_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Keep TIFF    : $([ -n "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')"
+echo "  Qualité image : $([ -n "$QUALITY" ] && echo "$QUALITY" || echo '98')"
+echo "  Format      : $([ -n "$FORMAT_FLAG" ] && echo "${FORMAT_FLAG#--format }" || echo 'avif')"
 echo "  Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')"
+if [ -n "$ONLY_FLAG" ]; then
+    echo "  Visualisations: uniquement${ONLY_FLAG#--only}"
+elif [ -n "$SKIP_FLAG" ]; then
+    echo "  Visualisations: tout sauf${SKIP_FLAG#--skip}"
+fi
 if [ -n "$FILE_ARGS" ]; then
     echo "  Fichiers   :${FILE_ARGS}"
 fi
 echo "============================================"
 
-CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG"
+CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG $QUALITY $FORMAT_FLAG"
 if [ -n "$GROUND_METHOD" ]; then
     CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD"
 fi
+if [ -n "$ONLY_FLAG" ]; then
+    CMD_ARGS="$CMD_ARGS $ONLY_FLAG"
+fi
+if [ -n "$SKIP_FLAG" ]; then
+    CMD_ARGS="$CMD_ARGS $SKIP_FLAG"
+fi
 if [ -n "$FILE_ARGS" ]; then
     CMD_ARGS="$CMD_ARGS --file $FILE_ARGS"
 fi