diff --git a/Dockerfile b/Dockerfile index 6ed79fc..d03c75b 100644 --- a/Dockerfile +++ b/Dockerfile @@ -14,7 +14,7 @@ RUN apt-get update && apt-get install -y --no-install-recommends \ wget \ && rm -rf /var/lib/apt/lists/* -WORKDIR /data +WORKDIR /app # Install Python packages via pip COPY requirements.txt . @@ -27,12 +27,18 @@ RUN pip3 install --no-cache-dir \ scikit-image \ scikit-learn \ scipy \ - tqdm + tqdm \ + Pillow # 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" -# Copy scripts +# Copy and install the pipeline package +COPY setup.py . +COPY lidar_pipeline/ ./lidar_pipeline/ +RUN pip3 install --no-cache-dir . + +# Copy backward-compatible entry point COPY process_lidar.py /usr/local/bin/ RUN chmod +x /usr/local/bin/process_lidar.py @@ -42,8 +48,10 @@ RUN groupadd -g 1000 lidar && \ mkdir -p /data/output /data/input && \ chown -R lidar:lidar /data /data/output /data/input +WORKDIR /data + USER lidar VOLUME ["/data"] -CMD ["python3", "/usr/local/bin/process_lidar.py", "/data/input", "-o", "/data/output"] \ No newline at end of file +CMD ["python3", "-m", "lidar_pipeline", "/data/input", "-o", "/data/output"] \ No newline at end of file diff --git a/README.md b/README.md index 6dd1db2..1b0de08 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ Workflow automatisé pour générer des visualisations exploitables à partir de données LiDAR pour la détection de structures archéologiques. -## Visualisations (21 par fichier) +## Visualisations (20 par fichier) ### Visualisations principales | # | Visualisation | Utilité archéologique | @@ -11,31 +11,47 @@ Workflow automatisé pour générer des visualisations exploitables à partir de | 2 | **Pente (Slope)** | Murs de soutènement, talus, changements brusques | | 3 | **Aspect (Orientation)** | Direction des pentes, exposition | | 4 | **Courbure (Curvature)** | Fossés, terrasses, talus, concavité/convexité | -| 5 | **Sky-View Factor** | Structures, tumulus, fondations (ray-tracing 16 azimuts) | -| 6 | **Local Relief Model** | Micro-reliefs, fossés, levées de terrain | -| 7 | **Positive Openness** | Élévations, tumulus, bâtiments (ray-tracing 8 directions) | -| 8 | **Negative Openness** | Cavités, fossés, souterrains (ray-tracing 8 directions) | +| 5 | **Éclairage solaire** | Simulation de l'éclairage matinal | +| 6 | **Sky-View Factor** | Structures, tumulus, fondations (ray-tracing 16 azimuts) | +| 7 | **Local Relief Model** | Micro-reliefs, fossés, levées de terrain | +| 8 | **Positive Openness** | Élévations, tumulus, bâtiments (ray-tracing 8 directions) | +| 9 | **Negative Openness** | Cavités, fossés, souterrains (ray-tracing 8 directions) | ### Visualisations avancées | # | Visualisation | Description | Détection | |---|--------------|-------------|-----------| -| 9 | **MSRM** | Multi-Scale Relief Model (sigma 5/10/25/50/100m) | Tumulus, fossés, murs à toutes échelles | -| 10 | **TPI multi-échelle** | Topographic Position Index (5m + 100m) | Crêtes, vallées, plateformes | -| 11 | **VAT composite** | Fusion hillshade+pente+SVF en RGB | Meilleure carte unique archéologique | +| 10 | **MSRM** | Multi-Scale Relief Model (sigma 5/10/25/50/100m) | Tumulus, fossés, murs à toutes les échelles | +| 11 | **TPI multi-échelle** | Topographic Position Index (5m + 100m) | Crêtes, vallées, plateformes | | 12 | **Dépressions** | Remplissage cuvettes + différence | Dolines, sinkholes, zones inondables | | 13 | **SAILORE** | LRM adaptatif (noyau = f(pente)) | Terrain hétérogène, tout relief | -| 14 | **Geomorphons** | 10 formes de terrain | Pics, crêtes, vallées, fosses, plateaux | -| 15 | **Rugosité** | Écart-type de l'élévation | Surfaces anthropiques vs naturelles | -| 16 | **Anomalies statistiques** | Z-score + Local Moran's I | Anomalies topographiques significatives | -| 17 | **Ondelette Mexican Hat** | CWT 2D multi-échelle | Tumulus, fossés circulaires | -| 18 | **Texture GLCM** | Contraste, entropie, homogénéité | Labour, surfaces archéologiques | -| 19 | **Accumulation de flux** | Algorithme D8 hydrologique | Fossés d'enceinte, routes antiques | +| 14 | **Rugosité** | Écart-type de l'élévation | Surfaces anthropiques vs naturelles | +| 15 | **Anomalies statistiques** | Z-score + Local Moran's I | Anomalies topographiques significatives | +| 16 | **Ondelette Mexican Hat** | CWT 2D multi-échelle | Tumulus, fossés circulaires | +| 17 | **Texture GLCM** | Contraste, entropie, homogénéité | Labour, surfaces archéologiques | +| 18 | **Accumulation de flux** | Algorithme D8 hydrologique | Fossés d'enceinte, routes antiques | ### Cartes de référence IGN | # | Visualisation | Source | |---|--------------|--------| -| 20 | **Photographie aérienne IGN** | Orthophotographie WMTS | -| 21 | **Carte topographique IGN** | Plan IGN V2 WMTS | +| 19 | **Photographie aérienne IGN** | Orthophotographie WMTS | +| 20 | **Carte topographique IGN** | Plan IGN V2 WMTS | + +## Architecture modulaire + +``` +lidar_pipeline/ +├── __init__.py # Exports publics +├── __main__.py # Point d'entrée: python -m lidar_pipeline +├── cli.py # argparse + logging + main() +├── gpu.py # CuPy/numpy abstraction (HAS_GPU, to_gpu, to_cpu, xp_*) +├── dtm.py # Classification PDAL + génération DTM +├── visualizations.py # Fonctions generate_* (20 visualisations) +├── ign.py # Téléchargement tuiles IGN + overlay +├── rendering.py # Colormaps, tif_to_png, rapport PDF +└── pipeline.py # LidarArchaeoPipeline (orchestration + registry) +``` + +Ajouter une visualisation = 1 fonction + 1 entrée dans `VIZ_STEPS` + 1 entrée dans `COLORMAPS`. ## Installation Docker @@ -46,86 +62,111 @@ mkdir -p input # Copiez vos fichiers .laz dans input/ cp /chemin/vos/fichiers/*.laz input/ -# Build l'image Docker (avec support GPU NVIDIA) -docker build -t lidar-archeo . +# Build l'image Docker +docker build -t lidar-lidar . ``` ## Utilisation -### Traitement avec accélération GPU (recommandé) +### Traitement complet avec GPU (recommandé) ```bash -# Nécessite une carte NVIDIA + nvidia-container-toolkit -docker run --rm --gpus all \ - -v $(pwd)/input:/data/input:ro \ - -v $(pwd)/output:/data/output \ - lidar-archeo +./run.sh -g ``` -### Traitement standard (CPU seul, sans GPU) +### Traitement standard (CPU seul) ```bash -docker run --rm \ - -v $(pwd)/input:/data/input:ro \ - -v $(pwd)/output:/data/output \ - lidar-archeo +./run.sh ``` -### Traitement parallèle multi-CPU -```bash -# Utiliser 4 CPU pour traiter plusieurs fichiers en parallèle -docker run --rm \ - -v $(pwd)/input:/data/input:ro \ - -v $(pwd)/output:/data/output \ - lidar-archeo \ - process_lidar.py /data/input -o /data/output -w 4 +### Options du script run.sh +``` +./run.sh [options] + -r RESOLUTION Résolution en m/px (défaut: 0.5) + -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 les WebP existent + --file NOM Traiter un seul fichier LAZ (pour tests) + -h Afficher l'aide ``` -### Résolution personnalisée +### Exemples ```bash -# Résolution fine (0.2m) - bâtiments individuels -docker run --rm \ - -v $(pwd)/input:/data/input:ro \ - -v $(pwd)/output:/data/output \ - lidar-archeo \ - process_lidar.py /data/input -o /data/output -r 0.2 +# Traitement standard +./run.sh -# Résolution standard (0.5m) - recommandée archéologie -docker run --rm \ - -v $(pwd)/input:/data/input:ro \ - -v $(pwd)/output:/data/output \ - lidar-archeo +# Avec accélération GPU +./run.sh -g + +# GPU + mode verbeux +./run.sh -g -v + +# GPU + 4 workers parallèles +./run.sh -g -w 4 + +# Haute résolution + GPU + debug +./run.sh -g -v -r 0.2 + +# Forcer la régénération de tous les fichiers +./run.sh -g --force + +# Traiter un seul fichier pour test +./run.sh -g --file 6881 ``` -### Combinaison résolution + multi-CPU +### Utilisation directe Docker ```bash -docker run --rm \ - --memory=16g \ - -v $(pwd)/input:/data/input:ro \ - -v $(pwd)/output:/data/output \ - lidar-archeo \ - process_lidar.py /data/input -o /data/output -r 0.5 -w 4 +# Traitement standard +docker run --rm -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data/output lidar-lidar + +# Avec GPU +docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data/output lidar-lidar + +# Mode verbeux +docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data/output \ + lidar-lidar python3 -m lidar_pipeline /data/input -o /data/output -v + +# Un seul fichier (test rapide) +docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data/output \ + lidar-lidar python3 -m lidar_pipeline /data/input -o /data/output --file 6881 + +# Forcer la régénération +docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data/output \ + lidar-lidar python3 -m lidar_pipeline /data/input -o /data/output --force ``` ## Structure des dossiers ``` . -├── input/ # Vos fichiers .laz (monté en read-only) +├── input/ # Fichiers .laz (monté en read-only) ├── output/ # Résultats générés │ ├── DTM/ # Modèles numériques de terrain (GeoTIFF) -│ ├── visualisations/ # Images PNG par fichier LAZ +│ ├── visualisations/ # Images WebP par fichier LAZ │ │ ├── fichier_6881/ # Un sous-dossier par fichier LAZ -│ │ │ ├── ..._hillshade_multi.png -│ │ │ ├── ..._svf.png -│ │ │ ├── ..._mslrm.png -│ │ │ └── ... (21 visualisations) +│ │ │ ├── ..._hillshade_multi.webp +│ │ │ ├── ..._svf.webp +│ │ │ ├── ..._mslrm.webp +│ │ │ └── ... (20 visualisations) │ │ └── fichier_6882/ │ │ └── ... │ └── rapports/ # Rapports PDF A3 par fichier │ ├── fichier_6881_rapport.pdf │ └── fichier_6882_rapport.pdf +├── lidar_pipeline/ # Package Python modulaire +│ ├── __init__.py +│ ├── cli.py +│ ├── gpu.py +│ ├── dtm.py +│ ├── visualizations.py +│ ├── ign.py +│ ├── rendering.py +│ └── pipeline.py +├── process_lidar.py # Point d'entrée compatible ├── Dockerfile -├── docker-compose.yml -└── process_lidar.py +├── run.sh +└── README.md ``` ## Paramètres @@ -135,36 +176,34 @@ docker run --rm \ | Résolution | `-r` | 0.5 | Résolution en mètres par pixel | | Workers | `-w` | 1 | Nombre de CPU pour traitement parallèle | | Output | `-o` | /data/output | Dossier de sortie | +| Force | `-f/--force` | off | Régénérer même si les WebP existent | +| File | `--file` | tous | Traiter un seul fichier (pour tests) | +| Verbose | `-v` | off | Mode verbeux (timestamps + niveaux) | +| Debug | `--debug` | off | Mode debug (détails internes) | ### Résolution recommandée -- `0.2` - Très fine, bâtiments individuels (lent) -- `0.5` - Recommandée archéologie (équilibre vitesse/détail) -- `1.0` - Rapide, grandes structures uniquement +- `0.2` — Très fine, bâtiments individuels (lent) +- `0.5` — Recommandée archéologie (équilibre vitesse/détail) +- `1.0` — Rapide, grandes structures uniquement ## Interprétation archéologique ### Pour détecter les cavités et souterrains -1. **Negative Openness** - Zones sombres = creux profonds -2. **Dépressions** - Carte spécifique des dolines et sinkholes -3. **Local Relief Model** - Zones bleues = dépressions -4. **Hillshade** - Ombres inhabituelles en forme de trous +1. **Negative Openness** — Zones sombres = creux profonds +2. **Dépressions** — Carte spécifique des dolines et sinkholes +3. **Local Relief Model** — Zones bleues = dépressions +4. **Hillshade** — Ombres inhabituelles en forme de trous ### Pour détecter structures et bâtiments anciens -1. **VAT composite** - Meilleure carte unique combinant hillshade+pente+SVF -2. **Sky-View Factor** - Structures géométriques claires -3. **MSRM** - Détection multi-échelle de tous les reliefs -4. **Geomorphons** - Classification automatique des formes - -### Pour anomalies statistiques -1. **Anomalies statistiques** - Zones significativement différentes du terrain -2. **Ondelette Mexican Hat** - Structures circulaires (tumulus, enclos) -3. **Rugosité** - Surfaces anthropiques vs naturelles -4. **Texture GLCM** - Labour ancien, chemins, murs enfouis +1. **MSRM** — Détection multi-échelle de tous les reliefs +2. **Sky-View Factor** — Structures géométriques claires +3. **SAILORE** — LRM adaptatif pour terrain hétérogène +4. **Anomalies statistiques** — Anomalies topographiques significatives ### Pour hydrologie et fossés -1. **Accumulation de flux** - Fossés d'enceinte, routes antiques -2. **Dépressions** - Zones de collecte d'eau, dolines -3. **Negative Openness** - Fossés et tranchées +1. **Accumulation de flux** — Fossés d'enceinte, routes antiques +2. **Dépressions** — Zones de collecte d'eau, dolines +3. **Negative Openness** — Fossés et tranchées ## Dépannage @@ -173,23 +212,15 @@ docker run --rm \ docker --version # Shell dans le conteneur -docker run --rm -it \ - -v $(pwd)/input:/data/input \ - -v $(pwd)/output:/data/output \ - --entrypoint bash lidar-archeo +docker run --rm -it -v $(pwd)/input:/data/input -v $(pwd)/output:/data/output \ + --entrypoint bash lidar-lidar # Reconstruire l'image -docker build --no-cache -t lidar-archeo . +docker build --no-cache -t lidar-lidar . # Nettoyer docker system prune -a ``` ### Erreur mémoire -Augmenter la mémoire Docker à 16Go+ pour les gros fichiers LiDAR HD. - -## Ressources - -- [PDAL Documentation](https://pdal.io/) -- [LiDAR Archéologie - Méthodes avancées](https://archaeologydataservice.ac.uk/) -- [Relief Visualization Toolbox](https://rvtpy.readthedocs.io/) \ No newline at end of file +Augmenter la mémoire Docker à 16Go+ pour les gros fichiers LiDAR HD. \ No newline at end of file diff --git a/lidar_pipeline/__init__.py b/lidar_pipeline/__init__.py new file mode 100644 index 0000000..599e80a --- /dev/null +++ b/lidar_pipeline/__init__.py @@ -0,0 +1,21 @@ +"""LiDAR Archaeological Pipeline — detection of archaeological structures from LiDAR data. + +Usage: + python -m lidar_pipeline /path/to/input -o /path/to/output -r 0.5 -v +""" + +# Lazy imports to avoid importing heavy dependencies at package import time + + +def __getattr__(name): + """Lazy-load heavy submodules only when accessed.""" + if name == 'LidarArchaeoPipeline': + from .pipeline import LidarArchaeoPipeline + return LidarArchaeoPipeline + if name == 'VIZ_STEPS': + from .pipeline import VIZ_STEPS + return VIZ_STEPS + if name == 'main': + from .cli import main + return main + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") \ No newline at end of file diff --git a/lidar_pipeline/__main__.py b/lidar_pipeline/__main__.py new file mode 100644 index 0000000..0bb4869 --- /dev/null +++ b/lidar_pipeline/__main__.py @@ -0,0 +1,6 @@ +"""Entry point for running the pipeline as: python -m lidar_pipeline""" + +from .cli import main + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py new file mode 100644 index 0000000..27e5263 --- /dev/null +++ b/lidar_pipeline/cli.py @@ -0,0 +1,155 @@ +"""Command-line interface for the LiDAR archaeological pipeline. + +Handles argument parsing, logging configuration, and entry point. +""" + +import argparse +import logging +import sys + +from .pipeline import LidarArchaeoPipeline +from .gpu import log_gpu_status + +logger = logging.getLogger("lidar") + + +def setup_logging(verbose=False, debug=False): + """Configure the 'lidar' logger. + + Args: + verbose: If True, include timestamps and level names. + debug: If True, set level to DEBUG and add file:line info. + """ + if debug: + level = logging.DEBUG + fmt = "%(asctime)s.%(msecs)03d %(levelname)-5s [%(filename)s:%(lineno)d] %(message)s" + elif verbose: + level = logging.INFO + fmt = "%(asctime)s %(levelname)-5s %(message)s" + else: + level = logging.INFO + fmt = "%(message)s" + + handler = logging.StreamHandler(sys.stdout) + handler.setFormatter(logging.Formatter(fmt, datefmt="%H:%M:%S")) + logger.setLevel(level) + logger.handlers.clear() + logger.addHandler(handler) + return logger + + +def main(): + """Entry point for the LiDAR archaeological pipeline.""" + parser = argparse.ArgumentParser( + description="Pipeline LiDAR pour détection archéologique", + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog="""\ +Exemples: + Traitement standard: + python -m lidar_pipeline /data/input -o /data/output + + Haute résolution avec accélération GPU: + python -m lidar_pipeline /data/input -o /data/output -r 0.2 -g + + Mode verbeux (timestamps): + python -m lidar_pipeline /data/input -o /data/output -v + + Mode debug (détails internes): + python -m lidar_pipeline /data/input -o /data/output --debug + + Forcer la régénération de tous les fichiers: + python -m lidar_pipeline /data/input -o /data/output --force + + Traiter un seul fichier (pour tests): + python -m lidar_pipeline /data/input -o /data/output --file LHD_FXX_1000_6881_PTS_LAMB93_IGN69.copc.laz + + Traitement parallèle (4 workers): + python -m lidar_pipeline /data/input -o /data/output -w 4 +""" + ) + parser.add_argument( + "input", + help="Dossier contenant les fichiers LAZ/LAS" + ) + parser.add_argument( + "-o", "--output", + default="/data/output", + help="Dossier de sortie (défaut: /data/output)" + ) + parser.add_argument( + "-r", "--resolution", + type=float, + default=0.5, + help="Résolution en mètres par pixel (défaut: 0.5)" + ) + parser.add_argument( + "-w", "--workers", + type=int, + default=1, + help="Nombre de workers pour traitement parallèle (défaut: 1)" + ) + parser.add_argument( + "-f", "--force", + action="store_true", + help="Régénérer tous les fichiers même si les WebP existent déjà" + ) + parser.add_argument( + "--file", + type=str, + default=None, + help="Traiter un seul fichier LAZ/LAS (pour tests, par nom partiel ou complet)" + ) + parser.add_argument( + "-v", "--verbose", + action="store_true", + help="Mode verbeux : affiche les timestamps et niveaux" + ) + parser.add_argument( + "--debug", + action="store_true", + help="Mode debug : affiche les détails internes (fichier:ligne)" + ) + + args = parser.parse_args() + + # Configure logging before any other output + setup_logging(verbose=args.verbose, debug=args.debug) + + logger.info("=" * 60) + logger.info("Pipeline LiDAR Archéologique") + logger.info("=" * 60) + + log_gpu_status() + + try: + pipeline = LidarArchaeoPipeline( + input_dir=args.input, + output_dir=args.output, + resolution=args.resolution, + workers=args.workers, + force=args.force + ) + + # If --file is specified, process only that single file + if args.file: + from pathlib import Path + input_dir = Path(args.input) + # Find matching file + matches = list(input_dir.glob(f"*{args.file}*")) + list(input_dir.glob(f"*{args.file}*.laz")) + list(input_dir.glob(f"*{args.file}*.las")) + # Remove duplicates + matches = list(dict.fromkeys(matches)) + if not matches: + logger.error(f"Aucun fichier trouvé pour: {args.file}") + sys.exit(1) + if len(matches) > 1: + logger.info(f"Plusieurs fichiers correspondent, utilisation du premier:") + for m in matches: + logger.info(f" {m.name}") + laz_file = matches[0] + logger.info(f"Traitement du fichier: {laz_file.name}") + pipeline.process_file(laz_file) + else: + pipeline.process_all() + except Exception as e: + logger.error(f"Erreur fatale: {e}", exc_info=True) + sys.exit(1) \ No newline at end of file diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py new file mode 100644 index 0000000..133dbdb --- /dev/null +++ b/lidar_pipeline/dtm.py @@ -0,0 +1,209 @@ +"""DTM generation from classified LiDAR point clouds. + +Handles ground classification via PDAL/SMRF and DTM rasterisation +using scipy binned_statistic_2d with gap-filling interpolation. +""" + +import json +import logging +import subprocess +from pathlib import Path + +import numpy as np +import rasterio +from rasterio.transform import from_bounds +from scipy import ndimage as scipy_ndimage +from scipy.ndimage import distance_transform_edt, gaussian_filter +from scipy.stats import binned_statistic_2d +from scipy import interpolate + +logger = logging.getLogger("lidar") + + +def create_smrf_pipeline(input_laz, output_las): + """Create a PDAL pipeline JSON for SMRF ground classification. + + Args: + input_laz: Path to input LAZ/LAS file. + output_las: Path to output classified LAS file. + + Returns: + JSON string of the PDAL pipeline. + """ + pipeline = { + "pipeline": [ + str(input_laz), + { + "type": "filters.smrf", + "ignore": "Classification[7:7]", + "slope": 1.0, + "window": 16.0, + "threshold": 0.5, + "scalar": 1.25 + }, + { + "type": "filters.range", + "limits": "Classification[2:2]" + }, + { + "type": "writers.las", + "filename": str(output_las), + "extra_dims": "all" + } + ] + } + return json.dumps(pipeline) + + +def classify_ground(laz_file, temp_dir): + """Classify ground points using PDAL SMRF filter. + + Args: + laz_file: Path to input LAZ/LAS file. + temp_dir: Directory for temporary files (pipeline.json, ground.las). + + Returns: + Path to classified ground LAS file, or None on failure. + """ + import laspy # noqa: ensure available + + output_las = temp_dir / f"{laz_file.stem}_ground.las" + + if output_las.exists(): + logger.info(" Classification déjà effectuée — fichier existant réutilisé") + return output_las + + pipeline_json = create_smrf_pipeline(laz_file, output_las) + pipeline_file = temp_dir / "pipeline.json" + + with open(pipeline_file, 'w') as f: + f.write(pipeline_json) + + try: + subprocess.run( + ["pdal", "pipeline", str(pipeline_file)], + capture_output=True, check=True + ) + logger.info(" ✓ Classification sol terminée") + return output_las + except subprocess.CalledProcessError as e: + error_msg = e.stderr.decode() + logger.error(f" ✗ Erreur PDAL: {error_msg}") + + # If error is about ReturnNumber=0, try filtering those points + if "ReturnNumber" in error_msg and "NumberOfReturns" in error_msg: + logger.info(" → Tentative de filtrage des points ReturnNumber=0...") + + filtered_pipeline = [ + {"type": "readers.las", "filename": str(laz_file)}, + {"type": "filters.range", "limits": "ReturnNumber[1:],NumberOfReturns[1:]"}, + {"type": "filters.smrf", "scalar": 1.25}, + {"type": "filters.range", "limits": "Classification[2:2]"}, + {"type": "writers.las", "filename": str(output_las), "extra_dims": "all"} + ] + + filtered_json = json.dumps(filtered_pipeline) + with open(pipeline_file, 'w') as f: + f.write(filtered_json) + + try: + subprocess.run( + ["pdal", "pipeline", str(pipeline_file)], + capture_output=True, check=True + ) + logger.info(" ✓ Classification sol terminée (points filtrés)") + return output_las + except subprocess.CalledProcessError as e2: + logger.error(f" ✗ Erreur même avec filtrage: {e2.stderr.decode()}") + return None + + return None + + +def create_dtm_fast(las_file, basename, dtm_dir, resolution): + """Create DTM using fast binning method with gap filling. + + Args: + las_file: Path to classified ground LAS file. + basename: Base name for output file. + dtm_dir: Directory for output DTM GeoTIFF. + resolution: Grid resolution in meters per pixel. + + Returns: + Path to output DTM GeoTIFF, or None on failure. + """ + import laspy + + logger.info(" → Génération DTM...") + + try: + las = laspy.read(str(las_file)) + + min_x, max_x = float(las.header.min[0]), float(las.header.max[0]) + min_y, max_y = float(las.header.min[1]), float(las.header.max[1]) + + width = int(np.ceil((max_x - min_x) / resolution)) + height = int(np.ceil((max_y - min_y) / resolution)) + + logger.debug(f" Bounds: X[{min_x:.1f}, {max_x:.1f}] Y[{min_y:.1f}, {max_y:.1f}]") + logger.debug(f" Grid: {width}x{height} pixels ({len(las.points):,} points)") + logger.info(f" Rasterisation {width}x{height} ({len(las.points):,} points)...") + + stat = binned_statistic_2d( + las.x, las.y, las.z, + statistic='mean', + bins=[width, height], + range=[[min_x, max_x], [min_y, max_y]] + ) + + dtm = stat.statistic.T + dtm = dtm[::-1, :] # Flip Y so north is at top + + # Fill gaps using interpolation + mask = np.isnan(dtm) + + if np.any(mask): + logger.info(" Remplissage des zones vides...") + + dist_to_nearest = distance_transform_edt(mask) + max_fill_distance = max(20, int(10 / resolution)) + + y_coords, x_coords = np.where(~mask) + z_values = dtm[~mask] + + if len(y_coords) > 100: + interp = interpolate.NearestNDInterpolator( + np.column_stack((y_coords, x_coords)), + z_values + ) + + y_missing, x_missing = np.where(mask) + dtm[y_missing, x_missing] = interp(y_missing, x_missing) + + dtm_smooth = gaussian_filter(dtm, sigma=2.0) + + fill_mask = mask & (dist_to_nearest <= max_fill_distance) + dtm[fill_mask] = dtm_smooth[fill_mask] + + far_mask = mask & (dist_to_nearest > max_fill_distance) + dtm[far_mask] = np.nan + + # Save as GeoTIFF + output_tif = dtm_dir / f"{basename}_dtm.tif" + transform = from_bounds(min_x, min_y, max_x, max_y, width, height) + + with rasterio.open( + output_tif, 'w', + driver='GTiff', height=height, width=width, + count=1, dtype='float32', + crs='EPSG:2154', transform=transform, + compress='lzw' + ) as dst: + dst.write(dtm.astype('float32'), 1) + + logger.info(f" ✓ DTM créé: {output_tif.name}") + return output_tif + + except Exception as e: + logger.error(f" ✗ Erreur DTM: {e}", exc_info=True) + return None \ No newline at end of file diff --git a/lidar_pipeline/gpu.py b/lidar_pipeline/gpu.py new file mode 100644 index 0000000..8bc4369 --- /dev/null +++ b/lidar_pipeline/gpu.py @@ -0,0 +1,73 @@ +"""GPU acceleration helpers for LiDAR pipeline. + +Provides CuPy/numpy abstraction layer. If CuPy is available and a CUDA GPU +is detected, array operations are accelerated on the GPU. Otherwise, all +operations fall back to numpy/scipy on CPU. +""" + +import logging +import numpy as np +from scipy import ndimage + +logger = logging.getLogger("lidar") + +# GPU detection - must happen at import time +HAS_GPU = False +_gpu_name = None +_gpu_mem_gb = 0 +_xp = np # Default: CPU + +try: + import cupy as cp + import cupyx.scipy.ndimage as cp_ndimage + + _gpu_info = cp.cuda.runtime.getDeviceProperties(0) + _gpu_name = _gpu_info['name'].decode() if isinstance(_gpu_info['name'], bytes) else str(_gpu_info['name']) + _gpu_mem_gb = _gpu_info['totalGlobalMem'] // (1024 ** 3) + HAS_GPU = True + _xp = cp +except (ImportError, Exception): + pass + + +def log_gpu_status(): + """Log GPU detection result. Called after logging is configured.""" + if HAS_GPU: + logger.info(f"GPU détectée: {_gpu_name} ({_gpu_mem_gb} Go VRAM)") + else: + logger.info("Pas de GPU — mode CPU uniquement") + + +def to_gpu(arr): + """Send array to GPU if available, otherwise return as float64 numpy.""" + if HAS_GPU: + return cp.asarray(arr.astype(np.float64)) + return arr.astype(np.float64) + + +def to_cpu(arr): + """Bring array back to CPU (numpy). No-op if already on CPU.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp.asnumpy(arr) + return arr + + +def xp_gaussian_filter(arr, sigma): + """Gaussian filter — uses GPU if array is on GPU, CPU otherwise.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp_ndimage.gaussian_filter(arr, sigma) + return ndimage.gaussian_filter(arr, sigma) + + +def xp_uniform_filter(arr, size): + """Uniform filter — uses GPU if array is on GPU, CPU otherwise.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp_ndimage.uniform_filter(arr, size) + return ndimage.uniform_filter(arr, size) + + +def xp_minimum_filter(arr, footprint=None, size=None): + """Minimum filter — uses GPU if array is on GPU, CPU otherwise.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp_ndimage.minimum_filter(arr, footprint=footprint, size=size) + return ndimage.minimum_filter(arr, footprint=footprint, size=size) \ No newline at end of file diff --git a/lidar_pipeline/ign.py b/lidar_pipeline/ign.py new file mode 100644 index 0000000..f401cc0 --- /dev/null +++ b/lidar_pipeline/ign.py @@ -0,0 +1,209 @@ +"""IGN (French National Geographic Institute) tile download and overlay generation. + +Downloads WMTS tiles from the IGN Geoportail and composites them into +GeoTIFF overlays matching the LiDAR DTM extent. +""" + +import logging +import math +import time +from pathlib import Path + +import numpy as np +import rasterio + +try: + from rasterio.warp import transform as warp_transform + HAS_WARP = True +except ImportError: + HAS_WARP = False + +logger = logging.getLogger("lidar") + + +def _lat_lon_to_tile(lat, lon, zoom): + """Convert lat/lon to Web Mercator tile coordinates.""" + n = 2 ** zoom + col = int((lon + 180) / 360 * n) + lat_rad = math.radians(lat) + row = int((1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n) + return col, row + + +def _lat_lon_to_px(lat, lon, zoom, tile_size=256): + """Convert lat/lon to Web Mercator pixel coordinates.""" + n = 2 ** zoom + px_x = (lon + 180) / 360 * n * tile_size + lat_rad = math.radians(lat) + px_y = (1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n * tile_size + return px_x, px_y + + +def download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom_level=15): + """Download IGN WMTS tiles for the given bounds using Web Mercator (PM). + + Args: + min_x, max_x, min_y, max_y: Bounds in Lambert 93. + layer: IGN WMTS layer name. + zoom_level: WMTS zoom level (default 15). + + Returns: + numpy array (H, W, 3) uint8, or None on failure. + """ + import urllib.request + import io + from PIL import Image as PILImage + + if not HAS_WARP: + logger.error(" ✗ rasterio.warp non disponible pour conversion de coordonnées") + return None + + try: + l93_xs = [min_x, max_x] + l93_ys = [max_y, min_y] + lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', l93_xs, l93_ys) + nw_lat, nw_lon = lats[0], lons[0] + se_lat, se_lon = lats[1], lons[1] + except Exception: + logger.error(" ✗ Impossible de convertir les coordonnées") + return None + + wmts_url = "https://data.geopf.fr/wmts" + tile_matrix_set = "PM" + tile_size = 256 + + col_min, row_min = _lat_lon_to_tile(nw_lat, nw_lon, zoom_level) + col_max, row_max = _lat_lon_to_tile(se_lat, se_lon, zoom_level) + + nw_px_x, nw_px_y = _lat_lon_to_px(nw_lat, nw_lon, zoom_level) + se_px_x, se_px_y = _lat_lon_to_px(se_lat, se_lon, zoom_level) + + out_width = int(se_px_x - nw_px_x) + out_height = int(se_px_y - nw_px_y) + + if out_width <= 0 or out_height <= 0 or out_width > 5000 or out_height > 5000: + return None + + composite = np.full((out_height, out_width, 3), 255, dtype=np.uint8) + + tiles_downloaded = 0 + fmt = "image/png" if 'PLAN' in layer else "image/jpeg" + + for col in range(col_min, col_max + 1): + for row in range(row_min, row_max + 1): + url = ( + f"{wmts_url}?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile" + f"&LAYER={layer}&STYLE=normal" + f"&TILEMATRIXSET={tile_matrix_set}" + f"&TILEMATRIX={zoom_level}&TILECOL={col}&TILEROW={row}" + f"&FORMAT={fmt}" + ) + + try: + req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'}) + with urllib.request.urlopen(req, timeout=10) as response: + tile_data = response.read() + tile_img = PILImage.open(io.BytesIO(tile_data)).convert('RGB') + tile_arr = np.array(tile_img) + + tile_origin_x = col * tile_size + tile_origin_y = row * tile_size + + px_x = int(tile_origin_x - nw_px_x) + px_y = int(tile_origin_y - nw_px_y) + + x_off = max(0, -px_x) + y_off = max(0, -px_y) + dst_x_start = max(0, px_x) + dst_y_start = max(0, px_y) + dst_x_end = min(out_width, px_x + tile_size) + dst_y_end = min(out_height, px_y + tile_size) + + src_x = x_off + src_y = y_off + src_w = dst_x_end - dst_x_start + src_h = dst_y_end - dst_y_start + + if src_w > 0 and src_h > 0 and tile_arr.shape[0] >= src_y + src_h and tile_arr.shape[1] >= src_x + src_w: + composite[dst_y_start:dst_y_end, dst_x_start:dst_x_end] = \ + tile_arr[src_y:src_y+src_h, src_x:src_x+src_w] + tiles_downloaded += 1 + + except Exception as e: + if tiles_downloaded == 0 and col == col_min and row == row_min: + logger.error(f" ✗ Erreur tuile IGN: {e}") + continue + + logger.info(f" → {tiles_downloaded} tuiles IGN téléchargées ({layer})") + if tiles_downloaded == 0: + return None + return composite + + +def generate_ign_overlay(dem_file, basename, vis_dir, resolution, layer, title, legend_label, description, out_suffix): + """Generate an IGN basemap overlay visualization. + + Args: + dem_file: Path to DTM GeoTIFF. + basename: Base name for output file. + vis_dir: Output directory for visualizations. + resolution: Grid resolution in m/px. + layer: IGN WMTS layer name. + title: Title for the overlay. + legend_label: Legend text. + description: Description text. + out_suffix: Suffix for output filename. + + Returns: + Path to output GeoTIFF, or None on failure. + """ + logger.info(f" → {title}...") + t0 = time.time() + + output = vis_dir / f"{basename}_{out_suffix}.tif" + + try: + with rasterio.open(dem_file) as src: + dem = src.read(1) + transform = src.transform + crs = src.crs + height, width = dem.shape + + top_left_x = transform.c + top_left_y = transform.f + pixel_size_x = transform.a + pixel_size_y = abs(transform.e) + + min_x = top_left_x + max_x = top_left_x + width * pixel_size_x + max_y = top_left_y + min_y = top_left_y - height * pixel_size_y + + zoom = 15 + + result = download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom) + + if result is None: + logger.error(" ✗ Impossible de télécharger les tuiles IGN") + return None + + from PIL import Image as PILImage + ign_pil = PILImage.fromarray(result) + ign_resized = ign_pil.resize((width, height), PILImage.LANCZOS) + ign_arr = np.array(ign_resized) + + with rasterio.open( + output, 'w', + driver='GTiff', height=height, width=width, + count=3, dtype='uint8', + crs=crs, transform=transform, + ) as dst: + for i in range(3): + dst.write(ign_arr[:, :, i], i + 1) + + logger.info(f" ✓ {title} terminé ({time.time()-t0:.1f}s)") + return output + + except Exception as e: + logger.error(f" ✗ Erreur IGN overlay: {e}", exc_info=True) + return None \ No newline at end of file diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py new file mode 100644 index 0000000..9fabb81 --- /dev/null +++ b/lidar_pipeline/pipeline.py @@ -0,0 +1,317 @@ +"""Pipeline orchestration for LiDAR archaeological analysis. + +LidarArchaeoPipeline coordinates the full processing chain: +1. Ground classification (PDAL/SMRF) +2. DTM generation +3. Visualization generation (19 products) +4. Rendering (WebP + PDF report) +""" + +import logging +import shutil +import time +from concurrent.futures import ProcessPoolExecutor, as_completed +from pathlib import Path +import subprocess + +from .dtm import classify_ground, create_dtm_fast +from .visualizations import ( + generate_hillshade, generate_slope, generate_aspect, generate_curvature, + generate_solar, generate_lrm, generate_svf, generate_openness, + generate_mslrm, generate_tpi, generate_depressions, generate_sailore, + generate_roughness, generate_anomalies, generate_wavelet, generate_texture, + generate_flow, +) +from .ign import generate_ign_overlay +from .rendering import tif_to_png, generate_pdf_report + +logger = logging.getLogger("lidar") + + +# Ordered list of visualization steps. +# Each entry: (name, function_or_lambda) +# Adding a new visualization = add a generate_* function + register here. +VIZ_STEPS = [ + ('hillshade', generate_hillshade), + ('slope', generate_slope), + ('aspect', generate_aspect), + ('curvature', generate_curvature), + ('solar', generate_solar), + ('svf', generate_svf), + ('lrm', generate_lrm), + ('pos_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=True)), + ('neg_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=False)), + ('mslrm', generate_mslrm), + ('tpi', generate_tpi), + ('depressions', generate_depressions), + ('sailore', generate_sailore), + ('roughness', generate_roughness), + ('anomalies', generate_anomalies), + ('wavelet', generate_wavelet), + ('texture', generate_texture), + ('flow', generate_flow), + ('ortho', lambda d, b, v, r: generate_ign_overlay( + d, b, v, r, + layer='ORTHOIMAGERY.ORTHOPHOTOS', + title='Photographie Aérienne IGN', + legend_label='Orthophotographie\nImage aérienne', + description='Photographie aérienne IGN (Orthophoto)', + out_suffix='ortho')), + ('topo', lambda d, b, v, r: generate_ign_overlay( + d, b, v, r, + layer='GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2', + title='Carte Topographique IGN', + legend_label='Carte IGN\nPlan topographique', + description='Carte topographique IGN (Plan IGN)', + out_suffix='topo')), +] + + +class LidarArchaeoPipeline: + """Orchestrates the LiDAR archaeological analysis pipeline.""" + + def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False): + self.input_dir = Path(input_dir) + self.output_dir = Path(output_dir) + self.resolution = resolution + self.workers = workers + self.force = force + self.temp_dir = self.output_dir / "temp" + + if not self.input_dir.exists(): + raise ValueError(f"Répertoire introuvable: {self.input_dir}") + + self.output_dir.mkdir(parents=True, exist_ok=True) + self.temp_dir.mkdir(exist_ok=True) + + 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]: + d.mkdir(exist_ok=True) + + 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") + logger.info(f" Workers : {workers}") + logger.info(f" Force : {'OUI' if self.force else 'non (skip existing)'}") + + def find_laz_files(self): + """Find all LAZ/LAS files in input directory.""" + files = list(self.input_dir.glob("*.laz")) + list(self.input_dir.glob("*.las")) + logger.info(f"{len(files)} fichier(s) LiDAR trouvé(s)") + for f in sorted(files): + logger.debug(f" {f.name}") + return sorted(files) + + def check_tools(self): + """Check that required external tools are available.""" + for name, cmd in [('pdal', 'pdal --version'), ('gdal', 'gdalinfo --version')]: + try: + result = subprocess.run(cmd.split(), capture_output=True, check=True, text=True) + version = result.stdout.strip().split('\n')[0] + logger.info(f" ✓ {name}: {version}") + except (subprocess.CalledProcessError, FileNotFoundError): + logger.error(f" ✗ {name} non disponible") + return False + return True + + def generate_all_visualizations(self, dtm_file, basename): + """Generate all archaeological visualizations for one DTM file. + + Returns a dict of {name: tif_path} for successful generations. + """ + logger.info(" Génération visualisations:") + + # Create per-file subdirectory + file_vis_dir = self.vis_dir / basename + file_vis_dir.mkdir(exist_ok=True) + + vis_results = {} + total = len(VIZ_STEPS) + + for idx, (name, func) in enumerate(VIZ_STEPS, 1): + # 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: + result = func(dtm_file, basename, file_vis_dir, self.resolution) + vis_results[name] = result + elapsed = time.time() - t0 + if result: + logger.info(f" [{idx}/{total}] ✓ {name} ({elapsed:.1f}s)") + else: + logger.warning(f" [{idx}/{total}] ✗ {name} — no output ({elapsed:.1f}s)") + except Exception as e: + vis_results[name] = None + logger.error(f" [{idx}/{total}] ✗ {name}: {e}", exc_info=True) + + # Convert to WebP (only newly generated TIFs, not skipped ones) + logger.info(" Conversion images WebP:") + 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, self.resolution) + if webp_file: + logger.info(f" ✓ {webp_file.name}") + + return vis_results + + def process_file(self, laz_file): + """Process a single LAZ file through the full pipeline.""" + basename = laz_file.stem + t_start = time.time() + + logger.info("=" * 60) + logger.info(f"FICHIER : {basename}") + logger.info("=" * 60) + + # Step 1: Ground classification + logger.info("[1/4] Classification du sol...") + t1 = time.time() + las_file = classify_ground(laz_file, self.temp_dir) + 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/4] Génération DTM...") + t2 = time.time() + dtm_file = create_dtm_fast(las_file, basename, self.dtm_dir, self.resolution) + 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)") + + # Step 3: Visualizations + logger.info("[3/4] Visualisations archéologiques...") + self.generate_all_visualizations(dtm_file, basename) + + # Step 4: PDF report + file_vis_dir = self.vis_dir / basename + logger.info("[4/4] Rapport PDF A3...") + t4 = time.time() + generate_pdf_report(basename, file_vis_dir, self.pdf_dir, self.resolution) + t_pdf = time.time() - t4 + logger.info(f" ✓ Rapport PDF terminé ({t_pdf:.1f}s)") + + 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") + return True + + def process_all(self): + """Process all LAZ files in input directory.""" + files = self.find_laz_files() + + if not files: + logger.error("Aucun fichier LAZ/LAS trouvé !") + return + + logger.info("=" * 60) + logger.info("PIPELINE ARCHÉOLOGIQUE LiDAR") + logger.info("=" * 60) + + logger.info("Vérification des outils...") + if not self.check_tools(): + logger.error("Outils manquants — abandon") + return + + results = {} + t_pipeline_start = time.time() + + if self.workers > 1 and len(files) > 1: + logger.info(f"Traitement parallèle avec {self.workers} workers...") + logger.info(f"Fichiers: {len(files)}") + with ProcessPoolExecutor(max_workers=self.workers) as executor: + future_to_file = { + executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution, self.force): laz_file + for laz_file in files + } + for idx, future in enumerate(as_completed(future_to_file), 1): + laz_file = future_to_file[future] + try: + success = future.result() + results[laz_file.name] = success + status = "✓" if success else "✗" + logger.info(f" [{idx}/{len(files)}] {status} {laz_file.name}") + except Exception as e: + logger.error(f" [{idx}/{len(files)}] ✗ {laz_file.name}: {e}") + logger.debug(f" Traceback:", exc_info=True) + results[laz_file.name] = False + 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 Exception as e: + logger.error(f"✗ Erreur traitement {laz_file.name}: {e}") + logger.debug("Traceback:", exc_info=True) + results[laz_file.name] = False + + # Summary + t_pipeline_total = time.time() - t_pipeline_start + success_count = sum(1 for v in results.values() if v) + fail_count = sum(1 for v in results.values() if not v) + + logger.info("=" * 60) + logger.info("RÉSUMÉ") + 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...") + try: + if self.temp_dir.exists(): + shutil.rmtree(self.temp_dir) + logger.info(" ✓ Fichiers temporaires supprimés") + except Exception as e: + 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): + """Standalone function for multiprocessing — creates its own pipeline instance. + + Each worker gets its own temp directory to avoid file conflicts. + """ + pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force) + basename = Path(laz_file_str).stem + pipeline.temp_dir = pipeline.output_dir / f"temp_{basename}" + pipeline.temp_dir.mkdir(exist_ok=True) + laz_file = Path(laz_file_str) + return pipeline.process_file(laz_file) \ No newline at end of file diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py new file mode 100644 index 0000000..9721b0c --- /dev/null +++ b/lidar_pipeline/rendering.py @@ -0,0 +1,605 @@ +"""Rendering module: colormap registry, GeoTIFF-to-WebP 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 +- generate_pdf_report(): generate an A3 PDF report with all visualizations +""" + +import logging +import time +from pathlib import Path + +import numpy as np +import rasterio + +try: + from rasterio.warp import transform as warp_transform + HAS_WARP = True +except ImportError: + HAS_WARP = False + +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rcParams +from matplotlib.gridspec import GridSpec +from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch +from mpl_toolkits.axes_grid1.inset_locator import inset_axes + +rcParams['figure.dpi'] = 150 +rcParams['savefig.dpi'] = 300 +rcParams['font.size'] = 10 + +logger = logging.getLogger("lidar") + + +# ============================================================ +# Colormap registry +# ============================================================ +# Each entry: keyword → (cmap, title, legend_label, description, vmin_mode, vmax_mode) +# vmin_mode/vmax_mode: 'percentile_X_Y' or '0_max_X' or 'symmetric_X_Y' or 'fixed_0_1' +# For RGB images (ortho/topo), special handling is done in tif_to_png. + +COLORMAPS = { + 'hillshade': { + 'cmap': 'gray', + 'title': 'Hillshade Multidirectionnel', + 'legend': 'Illumination combinée de 6 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre', + 'description': 'Ombres portées révélant micro-relief (murs, fossés, terrasses)', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'fixed', 'vmax_val': 1, + }, + 'slope': { + 'cmap': 'inferno', + 'title': 'Pente (Inclinaison du terrain)', + 'legend': 'Inclinaison en degrés\nMin: {vmin:.1f}° | Max: {vmax:.1f}°\nClair = Forte pente | Sombre = Terrain plat', + 'description': 'Murs, talus et bords ressortent en clair — terrain plat en sombre', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'percentile', 'vmax_pct': 95, + }, + 'aspect': { + 'cmap': 'hsv', + 'title': 'Aspect (Direction des pentes)', + 'legend': 'Direction vers laquelle le terrain descend\nRouge=Nord | Vert=Est | Cyan=Sud | Bleu=Ouest', + 'description': 'Orientation des pentes — utile pour distinguer structures des formes naturelles', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'fixed', 'vmax_val': 360, + }, + 'solar': { + 'cmap': 'gray', + 'title': 'Éclairage Solaire', + 'legend': 'Illumination solaire (azimut 90°, altitude 30°)\nClair = Face éclairée | Sombre = Zone d\'ombre', + 'description': 'Simulation de l\'éclairage solaire matinal', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'fixed', 'vmax_val': 1, + }, + 'curvature': { + 'cmap': 'RdYlBu_r', + 'title': 'Courbure (Convexité/Concavité du terrain)', + 'legend': 'Changement de pente (1/m)\nRouge = Convexe (sommet de mur, levée)\nBleu = Concave (fond de fossé, dépression)', + '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', + 'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'lrm': { + 'cmap': 'RdBu_r', + 'title': 'LRM - Local Relief Model (échelle unique 15m)', + 'legend': 'Écart local par rapport au terrain moyen (m)\nRouge = Surélévation (+{vmax:.2f}m)\nBleu = Dépression ({vmin:.2f}m)\nNoyau gaussien unique de 15m', + 'description': 'Micro-relief à 15m seulement — voir MSRM pour toutes les échelles', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'positive_openness': { + 'cmap': 'YlOrBr', + 'title': 'Openness Positive (ouverture vers le haut)', + 'legend': 'Angle d\'ouverture vers le haut (deg)\nClair = Vue dégagée vers le ciel (sommets, plateaux)\nSombre = Vue bloquée (vallées encaissées)', + 'description': 'Ray-tracing 8 directions — complémentaire de la négative pour détecter crêtes', + 'vmin_mode': 'percentile', 'vmin_pct': 10, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, + 'negative_openness': { + 'cmap': 'GnBu_r', + 'title': 'Openness Negative (ouverture vers le bas)', + 'legend': 'Angle d\'ouverture vers le bas (deg)\nClair = Surplomb (bords de fossé, grottes)\nSombre = Terrain plat (fonds de vallée)\nMeilleur détecteur de cavités et dolines', + 'description': 'Ray-tracing 8 directions — détecte fossés, dolines, souterrains', + 'vmin_mode': 'percentile', 'vmin_pct': 10, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, + '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', + 'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'depressions': { + 'cmap': 'YlOrRd', + 'title': 'Dépressions (Remplissage hydrologique)', + 'legend': 'Profondeur des cuvettes détectées (m)\nTransparent = Pas de dépression\nJaune = Dépression légère | Rouge = Dépression profonde\nMax: {vmax:.2f}m', + 'description': 'Simule le remplissage d\'eau — détecte dolines, sinkholes, cuvettes et zones inondables', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'percentile', 'vmax_pct': 99, + }, + 'sailore': { + 'cmap': 'seismic', + 'title': 'SAILORE - LRM Auto-Adaptatif', + 'legend': 'Relief local adaptatif (m)\nRouge = Surélévation | Bleu = Dépression\n\nDifférence avec LRM/MSRM:\nLRM = noyau fixe 15m\nMSRM = 5 noyaux fixes\nSAILORE = noyau adapté à la pente\nPlat=grand noyau | Pente=petit noyau', + 'description': 'Noyau qui s\'adapte à la pente locale — terrain plat=grand noyau, pente=petit noyau', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + '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', + '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)', + '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', + 'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'texture': { + 'cmap': 'inferno', + 'title': 'Texture GLCM (Contraste + Entropie - Homogénéité)', + 'legend': 'Analyse de la texture du relief (fenêtre 5m)\nClair = Texture hétérogène (labour, ruines, sol perturbé)\nSombre = Texture homogène (sol nu, route, zone plate)\n\nCombine contraste, entropie et homogénéité', + 'description': 'Distingue surfaces anthropiques (labour, chemins) des naturelles', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), + }, + 'flow': { + 'cmap': 'Blues', + 'title': 'Accumulation de Flux Hydrologique (D8)', + 'legend': 'Logarithme de l\'accumulation d\'eau\nBlanc = Pas de collecte (sommet, crête)\nBleu foncé = Collecte maximale (thalweg, fossé)\n\nSimule l\'écoulement de l\'eau de pluie\nDétecte fossés d\'enceinte, routes et drainage', + 'description': 'Algorithme D8 — simule le cheminement de l\'eau pour détecter fossés et routes antiques', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, +} + +# RGB entries (ortho/topo) are handled specially +RGB_LEGENDS = { + 'ortho': { + 'title': 'Photographie Aérienne IGN', + 'legend': 'Orthophotographie\nImage aérienne', + 'description': 'Photographie aérienne IGN (Orthophoto)', + }, + 'topo': { + 'title': 'Carte Topographique IGN', + 'legend': 'Carte IGN\nPlan topographique', + 'description': 'Carte topographique IGN (Plan IGN)', + }, +} + + +def _apply_colormap(data, tif_file): + """Apply the registered colormap normalization to data based on filename. + + Returns (data, cmap, title, legend_label, description, is_rgb). + """ + name = str(tif_file).lower() + + # Check for RGB first + for key in RGB_LEGENDS: + if key in name: + info = RGB_LEGENDS[key] + return data, None, info['title'], info['legend'], info['description'], True + + # Find matching colormap + for key, info in COLORMAPS.items(): + if key in name: + valid_data = data[~np.isnan(data)] if hasattr(data, 'compressed') else data.flatten() + valid_data = valid_data[~np.isnan(valid_data)] + + vmin = vmax = None + + # Compute vmin/vmax based on mode + if info['vmin_mode'] == 'fixed': + vmin = info['vmin_val'] + elif info['vmin_mode'] == 'percentile': + vmin = np.percentile(valid_data, info['vmin_pct']) + elif info['vmin_mode'] == 'symmetric': + vmax_abs = max(abs(np.percentile(valid_data, info['sym_pct'][0])), + abs(np.percentile(valid_data, info['sym_pct'][1])), 0.001) + vmin = -vmax_abs + vmax = vmax_abs # symmetric mode sets both vmin and vmax + + if vmax is None: + # Only compute vmax if not already set by symmetric mode + if info.get('vmax_mode') == 'fixed': + vmax = info['vmax_val'] + elif info.get('vmax_mode') == 'percentile': + vmax = np.percentile(valid_data, info['vmax_pct']) + elif info.get('vmax_mode') == 'symmetric': + vmax_abs = max(abs(np.percentile(valid_data, info['sym_pct'][0])), + abs(np.percentile(valid_data, info['sym_pct'][1])), 0.001) + vmax = vmax_abs + + # Apply normalization + if vmin is not None and vmax is not None: + data = np.clip((data - vmin) / max(vmax - vmin, 0.001), 0, 1) + + legend = info['legend'].format(vmin=vmin or 0, vmax=vmax or 0) + return data, info['cmap'], info['title'], legend, info['description'], False + + # Default: terrain colormap with percentile stretch + valid_data = data[~np.isnan(data)] if hasattr(data, 'compressed') else data.flatten() + valid_data = valid_data[~np.isnan(valid_data)] + p2, p98 = np.percentile(valid_data, (2, 98)) + data = np.clip((data - p2) / (p98 - p2), 0, 1) + title = Path(tif_file).stem.replace('_', ' ').title() + return data, 'terrain', title, 'Altitude normalisée', '', False + + +def tif_to_png(tif_file, vis_dir, resolution): + """Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar. + + Args: + tif_file: Path to input GeoTIFF. + vis_dir: Output directory for the WebP file. + resolution: Grid resolution in m/px. + + Returns: + Path to output WebP 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" + + try: + with rasterio.open(tif_file) as src: + is_rgb = src.count >= 3 and ('ortho' in str(tif_file) or 'topo' in str(tif_file)) + + if is_rgb: + data = src.read([1, 2, 3]) + data = np.moveaxis(data, 0, -1) + else: + data = src.read(1) + + nodata = src.nodata + transform = src.transform + crs = src.crs + + if is_rgb: + height, width, _ = data.shape + else: + height, width = data.shape + + top_left_x = transform.c + top_left_y = transform.f + pixel_size_x = transform.a + pixel_size_y = abs(transform.e) + + min_x = top_left_x + max_x = top_left_x + width * pixel_size_x + max_y = top_left_y + min_y = top_left_y - height * pixel_size_y + + # GPS coordinates + gps_coords = {} + if HAS_WARP and crs is not None: + try: + l93_xs = [min_x, max_x, min_x, max_x] + l93_ys = [max_y, max_y, min_y, min_y] + lons, lats = warp_transform(crs, 'EPSG:4326', l93_xs, l93_ys) + gps_coords = { + 'NW': (lats[0], lons[0]), + 'NE': (lats[1], lons[1]), + 'SW': (lats[2], lons[2]), + 'SE': (lats[3], lons[3]), + } + n_ticks = 5 + tick_l93_x = np.linspace(min_x, max_x, n_ticks) + tick_l93_y_bottom = np.full(n_ticks, min_y) + tick_lons, tick_lats = warp_transform(crs, 'EPSG:4326', tick_l93_x, tick_l93_y_bottom) + gps_coords['x_ticks'] = list(zip(tick_lons, tick_lats)) + tick_l93_y = np.linspace(min_y, max_y, n_ticks) + tick_l93_x_left = np.full(n_ticks, min_x) + tick_lons_y, tick_lats_y = warp_transform(crs, 'EPSG:4326', tick_l93_x_left, tick_l93_y) + gps_coords['y_ticks'] = list(zip(tick_lons_y, tick_lats_y)) + except Exception: + gps_coords = {} + + if nodata is not None and not is_rgb: + data = np.ma.masked_where((data == nodata) | np.isnan(data), data) + + if not is_rgb: + valid_data = data.compressed() if hasattr(data, 'compressed') else data.flatten() + valid_data = valid_data[~np.isnan(valid_data)] + + # Apply colormap + data, cmap, title, legend_label, description, is_rgb_result = _apply_colormap(data, tif_file) + + # Create figure + fig_width = 20 + map_aspect = height / width + fig = plt.figure(figsize=(fig_width, fig_width * map_aspect * 0.7 + 2.5), + facecolor='white') + gs = GridSpec(2, 1, height_ratios=[1.0, 0.06], + hspace=0.04, figure=fig, + left=0.06, right=0.88, top=0.93, bottom=0.08) + + ax = fig.add_subplot(gs[0]) + if is_rgb: + im = ax.imshow(data, aspect='equal', origin='upper') + else: + im = ax.imshow(data, cmap=cmap, aspect='equal', origin='upper') + + ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10) + + if not is_rgb: + cbar = plt.colorbar(im, ax=ax, pad=0.02, shrink=0.85, aspect=30) + cbar.ax.tick_params(labelsize=9, width=1.5) + cbar.outline.set_linewidth(1.5) + cbar.set_label(legend_label, fontsize=10, fontweight='bold') + else: + ax.text(1.02, 0.5, legend_label, transform=ax.transAxes, + fontsize=10, fontweight='bold', rotation=90, + verticalalignment='center', horizontalalignment='left') + + # GPS coordinate ticks + if gps_coords and 'x_ticks' in gps_coords: + x_pixel_positions = np.linspace(0, width - 1, len(gps_coords['x_ticks'])) + x_labels = [f"{lon:.5f}E" for lon, lat in gps_coords['x_ticks']] + ax.set_xticks(x_pixel_positions) + ax.set_xticklabels(x_labels, fontsize=7, rotation=30) + ax.set_xlabel('Longitude', fontsize=9, fontweight='bold') + + y_pixel_positions = np.linspace(0, height - 1, len(gps_coords['y_ticks'])) + y_labels = [f"{lat:.5f}N" for lon, lat in gps_coords['y_ticks']] + ax.set_yticks(y_pixel_positions) + ax.set_yticklabels(y_labels, fontsize=7) + ax.set_ylabel('Latitude', fontsize=9, fontweight='bold') + else: + x_ticks_count = 5 + x_positions = np.linspace(0, width - 1, x_ticks_count) + x_labels = [f"{(min_x + xp * pixel_size_x)/1000:.1f}" for xp in x_positions] + ax.set_xticks(x_positions) + ax.set_xticklabels(x_labels, fontsize=8) + ax.set_xlabel('Est (km) - Lambert 93', fontsize=9, fontweight='bold') + + y_ticks_count = 5 + y_positions = np.linspace(0, height - 1, y_ticks_count) + y_labels = [f"{(max_y - yp * pixel_size_y)/1000:.1f}" for yp in y_positions] + ax.set_yticks(y_positions) + ax.set_yticklabels(y_labels, fontsize=8) + ax.set_ylabel('Nord (km) - Lambert 93', fontsize=9, fontweight='bold') + + ax.tick_params(axis='both', which='both', direction='out', length=3, + width=0.8, colors='black') + for spine in ax.spines.values(): + spine.set_visible(True) + spine.set_color('black') + spine.set_linewidth(0.8) + + # North arrow + north_ax = inset_axes(ax, width="4%", height="7%", loc='upper right', + bbox_to_anchor=(-0.03, 0.12, 1, 1), bbox_transform=ax.transAxes) + north_ax.set_xlim(0, 1) + north_ax.set_ylim(0, 1) + north_ax.axis('off') + north_ax.plot([0.5, 0.5], [0.1, 0.65], color='black', linewidth=2.5, zorder=10) + north_ax.add_patch(MplPolygon([[0.5, 0.2], [0.35, 0.4], [0.5, 0.65], [0.65, 0.4]], + closed=True, facecolor='black', edgecolor='black', zorder=9)) + north_ax.text(0.5, 0.95, 'N', ha='center', va='top', + fontsize=13, fontweight='bold', color='black', zorder=11) + + # Bottom info bar + info_ax = fig.add_subplot(gs[1]) + info_ax.axis('off') + + extent_km_x = (max_x - min_x) / 1000 + extent_km_y = (max_y - min_y) / 1000 + + if is_rgb: + alt_min = alt_max = 0 + else: + alt_min = float(np.nanmin(valid_data)) if len(valid_data) > 0 else 0 + alt_max = float(np.nanmax(valid_data)) if len(valid_data) > 0 else 0 + + if gps_coords: + nw_lat, nw_lon = gps_coords['NW'] + se_lat, se_lon = gps_coords['SE'] + info_text = ( + f"GPS: {nw_lat:.5f}N {nw_lon:.5f}E - {se_lat:.5f}N {se_lon:.5f}E | " + f"EPSG:2154 | Res: {resolution}m/px | " + f"Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km" + ) + else: + info_text = ( + f"EPSG:2154 | X: {min_x:.0f}-{max_x:.0f} Y: {min_y:.0f}-{max_y:.0f} | " + f"Res: {resolution}m/px | Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km" + ) + + info_ax.text(0.01, 0.5, info_text, + transform=info_ax.transAxes, fontsize=8.5, + verticalalignment='center', family='monospace', + bbox=dict(boxstyle='round,pad=0.3', facecolor='#f0f0f0', + edgecolor='#aaaaaa', alpha=0.95)) + + # Scale bar + scale_m = 100 + pixels_per_meter = 1.0 / pixel_size_x + scale_px = int(scale_m * pixels_per_meter) + scale_bar_frac = scale_px / width + bar_left = 0.80 + bar_bottom = 0.15 + bar_width_frac = min(scale_bar_frac, 0.15) + bar_height = 0.35 + + info_ax.add_patch(RectPatch((bar_left, bar_bottom), bar_width_frac, bar_height, + facecolor='black', edgecolor='black', linewidth=0.5, + transform=info_ax.transAxes, clip_on=False)) + info_ax.text(bar_left + bar_width_frac / 2, bar_bottom + bar_height + 0.12, + f"{scale_m} m", ha='center', va='bottom', fontsize=9, fontweight='bold', + transform=info_ax.transAxes) + + fig.patch.set_facecolor('white') + + # Save as PNG then convert to WebP + png_temp = vis_dir / f"{tif_file.stem}_temp.png" + plt.savefig(png_temp, dpi=150, bbox_inches='tight', pad_inches=0.15, + facecolor='white', format='png') + plt.close() + + from PIL import Image as PILImage + img = PILImage.open(str(png_temp)) + img.save(str(webp_file), format='WEBP', lossless=True) + png_temp.unlink() + + # Delete source TIFF + tif_file.unlink() + + return webp_file + + except Exception as e: + logger.error(f" Erreur conversion WebP: {e}", exc_info=True) + return None + + +def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): + """Generate A3 PDF report for a LiDAR file with all visualizations. + + Page 1: Mise en situation (ortho + topo IGN side by side) + Pages 2+: Other visualizations (2 per page) + + Args: + basename: Base name for the report file. + vis_dir: Directory containing WebP visualization files. + pdf_dir: Directory for output PDF. + resolution: Grid resolution (used in info text). + + Returns: + Path to PDF file, or None on failure. + """ + from matplotlib.backends.backend_pdf import PdfPages + + pdf_file = pdf_dir / f"{basename}_rapport.pdf" + logger.info(f" → Génération rapport PDF A3: {pdf_file.name}") + t0 = time.time() + + # Look for WebPs in per-file subdirectory first, then fallback to main dir + file_vis_dir = vis_dir / basename + if file_vis_dir.exists(): + png_files = sorted(file_vis_dir.glob("*.webp")) + else: + png_files = sorted(vis_dir.glob(f"{basename}_*.webp")) + if not png_files: + logger.warning(f" ✗ Aucune image trouvée pour {basename}") + return None + + # Categorize + situ_files = [] + analysis_files = [] + + for f in png_files: + name = f.stem.lower() + if 'ortho' in name: + situ_files.insert(0, f) + elif 'topo' in name: + situ_files.append(f) + else: + analysis_files.append(f) + + # Sort analysis files by archaeological priority + order = ['mslrm', 'svf', 'negative_openness', + 'positive_openness', 'sailore', 'depressions', 'hillshade_multi', + 'lrm', 'tpi', 'slope', 'curvature', 'aspect', + 'roughness', 'anomalies', 'wavelet', 'texture', 'flow'] + + def sort_key(f): + name = f.stem.lower() + for i, key in enumerate(order): + if key in name: + return i + return len(order) + + analysis_files.sort(key=sort_key) + + a3_w, a3_h = 16.54, 11.69 + + try: + with PdfPages(str(pdf_file)) as pdf: + # Page 1: Mise en situation + if situ_files: + fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white') + n_situ = len(situ_files) + if n_situ == 2: + gs = fig.add_gridspec(1, 2, wspace=0.05, left=0.03, right=0.97, + top=0.92, bottom=0.06) + else: + gs = fig.add_gridspec(1, max(n_situ, 1), wspace=0.05, + left=0.03, right=0.97, top=0.92, bottom=0.06) + + fig.text(0.5, 0.97, f"Mise en situation - {basename}", + fontsize=20, fontweight='bold', ha='center', va='top') + + for i, f in enumerate(situ_files): + ax = fig.add_subplot(gs[0, i]) + img = plt.imread(str(f)) + ax.imshow(img) + ax.axis('off') + title = f.stem.replace(basename + '_', '').replace('_', ' ').title() + ax.set_title(title, fontsize=12, fontweight='bold', pad=5) + + pdf.savefig(fig, dpi=150) + plt.close(fig) + + # Pages 2+: Analysis maps (2 per page) + for page_start in range(0, len(analysis_files), 2): + page_files = analysis_files[page_start:page_start + 2] + + fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white') + + if len(page_files) == 2: + gs = fig.add_gridspec(1, 2, wspace=0.08, left=0.03, right=0.97, + top=0.93, bottom=0.05) + else: + gs = fig.add_gridspec(1, 1, left=0.05, right=0.95, + top=0.93, bottom=0.05) + + for i, f in enumerate(page_files): + ax = fig.add_subplot(gs[0, i]) + img = plt.imread(str(f)) + ax.imshow(img) + ax.axis('off') + title = f.stem.replace(basename + '_', '').replace('_', ' ').title() + ax.set_title(title, fontsize=11, fontweight='bold', pad=3) + + page_num = (page_start // 2) + 2 + fig.text(0.99, 0.01, f"Page {page_num}", fontsize=8, + ha='right', va='bottom', color='gray') + + pdf.savefig(fig, dpi=150) + plt.close(fig) + + logger.info(f" ✓ Rapport PDF terminé ({time.time()-t0:.1f}s)") + return pdf_file + + except Exception as e: + logger.error(f" ✗ Erreur PDF: {e}", exc_info=True) + return None \ No newline at end of file diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py new file mode 100644 index 0000000..b6dc4a4 --- /dev/null +++ b/lidar_pipeline/visualizations.py @@ -0,0 +1,740 @@ +"""Terrain visualization functions for LiDAR archaeological analysis. + +Each function takes (dem_file, basename, vis_dir, resolution) as explicit +parameters and returns the path to the output GeoTIFF file, or None on error. +""" + +import logging +import time +from pathlib import Path + +import numpy as np +import rasterio +from scipy import ndimage +from scipy.ndimage import generic_filter, gaussian_filter, uniform_filter, minimum_filter +from scipy.stats import binned_statistic_2d + +from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter + +logger = logging.getLogger("lidar") + +# Use CuPy array module when available +if HAS_GPU: + import cupy as cp + xp = cp +else: + xp = np + + +def _save_tif(output_path, data, transform, crs, dtype='float32', count=1): + """Helper to save a 2D or 3D array as GeoTIFF.""" + if data.ndim == 2: + height, width = data.shape + with rasterio.open( + output_path, 'w', driver='GTiff', + height=height, width=width, count=count, + dtype=dtype, crs=crs, transform=transform, compress='lzw' + ) as dst: + dst.write(data.astype(dtype), 1) + elif data.ndim == 3: + bands, height, width = data.shape + with rasterio.open( + output_path, 'w', driver='GTiff', + height=height, width=width, count=bands, + dtype=dtype, crs=crs, transform=transform, compress='lzw' + ) as dst: + for i in range(bands): + dst.write(data[i].astype(dtype), i + 1) + + +def _read_dem(dem_file): + """Read DEM file and return (data, transform, crs).""" + with rasterio.open(dem_file) as src: + return src.read(1), src.transform, src.crs + + +# ============================================================ +# Core terrain visualizations +# ============================================================ + +def generate_hillshade(dem_file, basename, vis_dir, resolution): + """Generate multi-directional hillshade (NW, NE, SW, SE).""" + logger.info(" → Hillshade multidirectionnel...") + t0 = time.time() + output = vis_dir / f"{basename}_hillshade_multi.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + dy, dx = np.gradient(dem) + + azimuts = [315, 45, 225, 135] + altitude = 30 + hillshades = [] + + for az in azimuts: + az_rad = np.radians(az) + alt_rad = np.radians(altitude) + slope = np.arctan(np.sqrt(dx**2 + dy**2)) + aspect = np.arctan2(dy, dx) + hs = np.sin(alt_rad) * np.sin(slope) + \ + np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect) + hillshades.append(np.clip(hs, 0, 1)) + + combined = np.mean(hillshades, axis=0) + _save_tif(output, combined, transform, crs) + logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur hillshade: {e}", exc_info=True) + return None + + +def generate_slope(dem_file, basename, vis_dir, resolution): + """Generate slope map (degrees).""" + logger.info(" → Pente (Slope)...") + t0 = time.time() + output = vis_dir / f"{basename}_slope.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + dy, dx = np.gradient(dem) + slope = np.arctan(np.sqrt(dx**2 + dy**2)) * 180 / np.pi + _save_tif(output, slope, transform, crs) + logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur slope: {e}", exc_info=True) + return None + + +def generate_aspect(dem_file, basename, vis_dir, resolution): + """Generate aspect (slope orientation) map.""" + logger.info(" → Aspect (Orientation)...") + t0 = time.time() + output = vis_dir / f"{basename}_aspect.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + dy, dx = np.gradient(dem) + aspect = np.arctan2(dy, dx) * 180 / np.pi + aspect = np.mod(aspect, 360) + _save_tif(output, aspect, transform, crs) + logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur aspect: {e}", exc_info=True) + return None + + +def generate_curvature(dem_file, basename, vis_dir, resolution): + """Generate curvature (terrain concavity/convexity) map.""" + logger.info(" → Courbure (Curvature)...") + t0 = time.time() + output = vis_dir / f"{basename}_curvature.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + dz_dx = np.gradient(dem, axis=1) + dz_dy = np.gradient(dem, axis=0) + d2z_dx2 = np.gradient(dz_dx, axis=1) + d2z_dy2 = np.gradient(dz_dy, axis=0) + curvature = (d2z_dx2 + d2z_dy2) / 2 + _save_tif(output, curvature, transform, crs) + logger.info(f" ✓ Courbure terminée ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur curvature: {e}", exc_info=True) + return None + + +def generate_solar(dem_file, basename, vis_dir, resolution): + """Generate solar irradiance simulation.""" + logger.info(" → Éclairage Solaire...") + t0 = time.time() + output = vis_dir / f"{basename}_solar.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + dy, dx = np.gradient(dem) + slope = np.arctan(np.sqrt(dx**2 + dy**2)) + aspect = np.arctan2(dy, dx) + az_rad = np.radians(90) + alt_rad = np.radians(30) + solar = np.sin(alt_rad) * np.sin(slope) + \ + np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect) + solar = np.clip(solar, 0, 1) + _save_tif(output, solar, transform, crs) + logger.info(f" ✓ Solaire terminé ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur solar: {e}", exc_info=True) + return None + + +# ============================================================ +# GPU-accelerated visualizations +# ============================================================ + +def generate_lrm(dem_file, basename, vis_dir, resolution): + """Local Relief Model - deviation from local mean (GPU if available).""" + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → Local Relief Model{gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_lrm.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + local_mean = xp_gaussian_filter(dem, sigma=15/resolution) + lrm = dem - local_mean + lrm_np = to_cpu(lrm).astype(np.float32) + _save_tif(output, lrm_np, transform, crs) + logger.info(f" ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur LRM: {e}", exc_info=True) + return None + + +def generate_svf(dem_file, basename, vis_dir, resolution): + """Sky-View Factor - ray-tracing on 16 azimuths (GPU if available). + + For each pixel, trace rays in N directions, find the max horizon + angle in each direction, then SVF = (1/N) * sum(cos²(horizon_angle)). + Valleys/crevices have low SVF (obstructed sky), ridges/peaks have high SVF. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → Sky-View Factor (ray-tracing){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_svf.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + rows, cols = dem_np.shape + res = resolution + + dem = to_gpu(dem_np) + + n_dirs = 16 + angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) + dx = np.cos(angles) + dy = np.sin(angles) + max_dist = int(50 / res) + + padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) + svf = xp.zeros_like(dem) + + for d_idx in range(n_dirs): + ddx, ddy = dx[d_idx], dy[d_idx] + horizon = xp.zeros_like(dem) + + for step in range(1, max_dist + 1): + px = int(round(ddx * step)) + py = int(round(ddy * step)) + dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) + if dist_m < res * 0.5: + continue + + elev_diff = padded[max_dist + py:max_dist + py + rows, + max_dist + px:max_dist + px + cols] - dem + angle = xp.arctan2(elev_diff, dist_m) + horizon = xp.where(xp.isnan(angle), horizon, + xp.maximum(horizon, xp.nan_to_num(angle, nan=0))) + + svf += xp.cos(xp.pi / 2 - horizon) ** 2 + + svf /= n_dirs + svf_np = to_cpu(svf).astype(np.float32) + _save_tif(output, svf_np, transform, crs) + logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur SVF: {e}", exc_info=True) + return None + + +def generate_openness(dem_file, basename, vis_dir, resolution, positive=True): + """Positive/Negative Openness - true zenith/nadir angle computation (GPU if available). + + For each pixel, in 8 directions (N, NE, E, SE, S, SW, W, NW): + - Positive openness: max zenith angle (angle from vertical to highest visible terrain) + - Negative openness: max nadir angle (angle from vertical down to lowest terrain) + Result is averaged across all 8 directions. + """ + name = "positive_openness" if positive else "negative_openness" + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_{name}.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + rows, cols = dem_np.shape + res = resolution + + dem = to_gpu(dem_np) + + n_dirs = 8 + angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) + dx = np.cos(angles) + dy = np.sin(angles) + max_dist = int(50 / res) + + padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) + openness_sum = xp.zeros_like(dem) + + for d_idx in range(n_dirs): + ddx, ddy = dx[d_idx], dy[d_idx] + max_angle = xp.zeros_like(dem) + + for step in range(1, max_dist + 1): + px = int(round(ddx * step)) + py = int(round(ddy * step)) + dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) + if dist_m < res * 0.5: + continue + + elev_diff = padded[max_dist + py:max_dist + py + rows, + max_dist + px:max_dist + px + cols] - dem + + if positive: + angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) + else: + angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) + + max_angle = xp.where(xp.isnan(angle), max_angle, + xp.maximum(max_angle, xp.nan_to_num(angle, nan=0))) + + openness_sum += max_angle + + openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32) + _save_tif(output, openness_result, transform, crs) + logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur openness: {e}", exc_info=True) + return None + + +def generate_mslrm(dem_file, basename, vis_dir, resolution): + """Multi-Scale Relief Model (MSRM) - LRM at 5 scales combined (GPU if available).""" + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_mslrm.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + + sigmas = [5, 10, 25, 50, 100] + lrm_stack = [] + + for sigma in sigmas: + sigma_px = sigma / resolution + local_mean = xp_gaussian_filter(dem, sigma=sigma_px) + lrm = dem - local_mean + lrm_norm = lrm / max(float(xp.nanstd(lrm)), 0.01) + lrm_stack.append(lrm_norm) + + mslrm = xp.sqrt(xp.mean(xp.array(lrm_stack) ** 2, axis=0)) + mslrm_np = to_cpu(mslrm).astype(np.float32) + _save_tif(output, mslrm_np, transform, crs) + logger.info(f" ✓ MSRM terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur MSRM: {e}", exc_info=True) + return None + + +def generate_tpi(dem_file, basename, vis_dir, resolution): + """Multi-Scale Topographic Position Index (GPU if available). + + TPI = elevation - mean(neighborhood). + Computed at fine (5m) and broad (100m) scales. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → TPI multi-échelle{gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_tpi.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + + fine_size = int(5 / resolution) + if fine_size % 2 == 0: + fine_size += 1 + tpi_fine = dem - xp_uniform_filter(dem, size=fine_size) + + broad_size = int(100 / resolution) + if broad_size % 2 == 0: + broad_size += 1 + tpi_broad = dem - xp_uniform_filter(dem, size=broad_size) + + fine_std = max(float(xp.nanstd(tpi_fine)), 0.01) + broad_std = max(float(xp.nanstd(tpi_broad)), 0.01) + tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std) + tpi_np = to_cpu(tpi_combined).astype(np.float32) + + _save_tif(output, tpi_np, transform, crs) + logger.info(f" ✓ TPI terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur TPI: {e}", exc_info=True) + return None + + +# ============================================================ +# Depression / hydrology +# ============================================================ + +def generate_depressions(dem_file, basename, vis_dir, resolution): + """Depression detection using hydrological sink filling.""" + logger.info(" → Détection dépressions (hydrologique)...") + t0 = time.time() + output = vis_dir / f"{basename}_depressions.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + + from scipy.ndimage import binary_dilation, generate_binary_structure + + dem_filled = dem.copy() + nodata_mask = np.isnan(dem_filled) + dem_filled[nodata_mask] = np.nanmax(dem) + 1000 + + struct = generate_binary_structure(2, 2) + changed = True + iterations = 0 + max_iter = 100 + + while changed and iterations < max_iter: + from scipy.ndimage import minimum_filter as scipy_min_filter + neighbor_min = scipy_min_filter(dem_filled, footprint=struct) + sinks = (dem_filled < neighbor_min) & ~nodata_mask + + if not np.any(sinks): + break + + new_dem = np.maximum(dem_filled, neighbor_min) + new_dem[nodata_mask] = np.nan + changed = np.any(new_dem != dem_filled) + dem_filled = new_dem + iterations += 1 + + depressions = dem_filled - dem + depressions[nodata_mask] = np.nan + depressions = np.where(depressions > 0.01, depressions, 0) + + _save_tif(output, depressions, transform, crs) + logger.info(f" ✓ Dépressions terminé ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur dépressions: {e}", exc_info=True) + return None + + +# ============================================================ +# SAILORE +# ============================================================ + +def generate_sailore(dem_file, basename, vis_dir, resolution): + """SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available). + + Kernel size adapts to local slope: flat areas get larger kernels, + steep areas get smaller kernels. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → SAILORE (LRM adaptatif){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_sailore.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + + gy, gx = xp.gradient(dem, resolution) + slope = xp.arctan(xp.sqrt(gx**2 + gy**2)) + slope_deg = xp.degrees(slope) + + sigma_min = 2.0 / resolution + sigma_max = 25.0 / resolution + slope_norm = xp.clip(slope_deg / 30.0, 0, 1) + adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min) + + lrm_fine = dem - xp_gaussian_filter(dem, sigma=sigma_min) + lrm_medium = dem - xp_gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2) + lrm_coarse = dem - xp_gaussian_filter(dem, sigma=sigma_max) + + w_fine = slope_norm + w_medium = 1 - 2 * xp.abs(slope_norm - 0.5) + w_coarse = 1 - slope_norm + w_total = w_fine + w_medium + w_coarse + w_total[w_total == 0] = 1 + + sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total + sailore_np = to_cpu(sailore).astype(np.float32) + + _save_tif(output, sailore_np, transform, crs) + logger.info(f" ✓ SAILORE terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur SAILORE: {e}", exc_info=True) + return None + + +# ============================================================ +# Roughness +# ============================================================ + +def generate_roughness(dem_file, basename, vis_dir, resolution): + """Surface roughness - standard deviation of elevation in a window.""" + logger.info(" → Rugosité de surface...") + t0 = time.time() + output = vis_dir / f"{basename}_roughness.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + window_size = int(5 / resolution) + if window_size % 2 == 0: + window_size += 1 + + def std_filter(arr): + return np.nanstd(arr) + + roughness = generic_filter(dem.astype(np.float64), std_filter, + size=window_size, mode='nearest') + _save_tif(output, roughness, transform, crs) + logger.info(f" ✓ Rugosité terminée ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur rugosité: {e}", exc_info=True) + return None + + +# ============================================================ +# Anomalies +# ============================================================ + +def generate_anomalies(dem_file, basename, vis_dir, resolution): + """Statistical anomaly detection - z-score of local relief + Local Moran's I.""" + logger.info(" → Détection anomalies statistiques...") + t0 = time.time() + output = vis_dir / f"{basename}_anomalies.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + + lrm = dem - gaussian_filter(dem, sigma=15 / resolution) + lrm_mean = np.nanmean(lrm) + lrm_std = max(np.nanstd(lrm), 0.01) + z_score = (lrm - lrm_mean) / lrm_std + + window = int(10 / resolution) + if window % 2 == 0: + window += 1 + + local_mean = uniform_filter(z_score, size=window) + morans_i = z_score * (local_mean - np.nanmean(z_score)) / np.nanstd(z_score) + anomaly_score = np.abs(z_score) * np.sign(morans_i) + + _save_tif(output, anomaly_score, transform, crs) + logger.info(f" ✓ Anomalies terminé ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur anomalies: {e}", exc_info=True) + return None + + +# ============================================================ +# Wavelet +# ============================================================ + +def generate_wavelet(dem_file, basename, vis_dir, resolution): + """Mexican Hat wavelet multi-scale analysis (GPU if available). + + CWT 2D at multiple scales to detect circular features. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → Ondelette Mexican Hat multi-échelle{gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_wavelet.tif" + + try: + dem_np, transform, crs = _read_dem(dem_file) + dem = to_gpu(dem_np) + + scales = [2, 5, 10, 20, 50] + wavelet_stack = [] + + for scale_m in scales: + sigma_px = scale_m / resolution + if HAS_GPU: + from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace + response = -gpu_gaussian_laplace(dem, sigma=sigma_px) + else: + from scipy.ndimage import gaussian_laplace + response = to_gpu(-gaussian_laplace(to_cpu(dem).astype(np.float64), sigma=sigma_px)) + response /= max(float(xp.nanstd(response)), 0.01) + wavelet_stack.append(response) + + combined = xp.sqrt(xp.mean(xp.array(wavelet_stack) ** 2, axis=0)) + combined_np = to_cpu(combined).astype(np.float32) + + _save_tif(output, combined_np, transform, crs) + logger.info(f" ✓ Ondelette terminée ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur ondelette: {e}", exc_info=True) + return None + + +# ============================================================ +# Texture GLCM +# ============================================================ + +def generate_texture(dem_file, basename, vis_dir, resolution): + """GLCM texture analysis on hillshade - contrast, entropy, homogeneity.""" + logger.info(" → Texture GLCM...") + t0 = time.time() + output = vis_dir / f"{basename}_texture.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + + gy, gx = np.gradient(dem, resolution) + slope = np.arctan(np.sqrt(gx**2 + gy**2)) + alt_rad = np.radians(45) + az_rad = np.radians(315) + shading = (np.sin(alt_rad) * np.cos(slope) + + np.cos(alt_rad) * np.sin(slope) * + np.cos(az_rad - np.arctan2(gy, gx))) + hillshade = np.clip(shading, 0, 1) + + valid = hillshade[~np.isnan(hillshade)] + if len(valid) == 0: + raise ValueError("No valid data for texture analysis") + lo, hi = np.percentile(valid, (1, 99)) + img = np.clip((hillshade - lo) / max(hi - lo, 0.001), 0, 1) + img_uint8 = (img * 255).astype(np.uint8) + + window = int(5 / resolution) + if window % 2 == 0: + window += 1 + + def local_variance(arr): + return np.var(arr.astype(np.float64)) + + def local_entropy(arr): + hist, _ = np.histogram(arr.astype(np.float64), bins=16, range=(0, 256)) + hist = hist / max(hist.sum(), 1) + hist = hist[hist > 0] + return -np.sum(hist * np.log2(hist)) + + def local_homogeneity(arr): + arr_f = arr.astype(np.float64) + return np.mean(1.0 / (1.0 + (arr_f - np.mean(arr_f)) ** 2)) + + contrast = generic_filter(img_uint8.astype(np.float64), local_variance, + size=window, mode='nearest') + entropy = generic_filter(img_uint8.astype(np.float64), local_entropy, + size=window, mode='nearest') + homogeneity = generic_filter(img_uint8.astype(np.float64), local_homogeneity, + size=window, mode='nearest') + + def norm(arr): + valid_arr = arr[~np.isnan(arr)] + if len(valid_arr) == 0: + return arr + std_val = max(np.std(valid_arr), 0.01) + return (arr - np.mean(valid_arr)) / std_val + + texture_combined = (0.4 * norm(contrast) + + 0.4 * norm(entropy) - + 0.2 * norm(homogeneity)) + + _save_tif(output, texture_combined, transform, crs) + logger.info(f" ✓ Texture terminée ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur texture GLCM: {e}", exc_info=True) + return None + + +# ============================================================ +# Flow accumulation +# ============================================================ + +def generate_flow(dem_file, basename, vis_dir, resolution): + """Flow accumulation using D8 algorithm. + + Identifies drainage patterns, ditches, and enclosure boundaries. + """ + logger.info(" → Accumulation de flux D8...") + t0 = time.time() + output = vis_dir / f"{basename}_flow.tif" + + try: + dem, transform, crs = _read_dem(dem_file) + rows, cols = dem.shape + nodata_mask = np.isnan(dem) + + from scipy.ndimage import minimum_filter as scipy_min_filter, generate_binary_structure + + dem_filled = dem.copy() + dem_filled[nodata_mask] = np.nanmax(dem) + 1000 + + struct = generate_binary_structure(2, 2) + for _ in range(50): + neighbor_min = scipy_min_filter(dem_filled, footprint=struct) + sinks = (dem_filled < neighbor_min) & ~nodata_mask + if not np.any(sinks): + break + dem_filled = np.where(sinks, neighbor_min, dem_filled) + + dem_filled[nodata_mask] = np.nan + + dx8 = [1, 1, 0, -1, -1, -1, 0, 1] + dy8 = [0, 1, 1, 1, 0, -1, -1, -1] + dist8 = [1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2)] + + flow_dir = np.full((rows, cols), -1, dtype=np.int8) + max_slope = np.full((rows, cols), 0.0) + + padded = np.pad(dem_filled, 1, mode='constant', + constant_values=np.nanmax(dem_filled) + 10000) + + for d in range(8): + nx = 1 + dx8[d] + ny = 1 + dy8[d] + neighbor_elev = padded[ny:ny + rows, nx:nx + cols] + slope = (dem_filled - neighbor_elev) / (dist8[d] * resolution) + slope[nodata_mask] = -1 + better = slope > max_slope + flow_dir[better] = d + max_slope[better] = slope[better] + + flat_dem = dem_filled[~nodata_mask].flatten() + valid_indices = np.where(~nodata_mask.flatten())[0] + sort_order = valid_indices[np.argsort(-flat_dem)] + + flow_acc = np.ones((rows, cols), dtype=np.float32) + flow_acc[nodata_mask] = 0 + + for idx in sort_order: + r, c = divmod(idx, cols) + d = flow_dir[r, c] + if d < 0: + continue + nr, nc = r + dy8[d], c + dx8[d] + if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]: + flow_acc[nr, nc] += flow_acc[r, c] + + flow_log = np.log1p(flow_acc) + _save_tif(output, flow_log, transform, crs) + logger.info(f" ✓ Flux terminé ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur flux: {e}", exc_info=True) + return None \ No newline at end of file diff --git a/process_lidar.py b/process_lidar.py index 8da10f7..5946283 100755 --- a/process_lidar.py +++ b/process_lidar.py @@ -1,2745 +1,13 @@ #!/usr/bin/env python3 -""" -Pipeline LiDAR pour détection archéologique -Visualisations générées avec numpy/cupy + rasterio/matplotlib -Support GPU via CuPy si disponible, fallback numpy sinon +"""Backward-compatible entry point for the LiDAR archaeological pipeline. + +This file exists for compatibility with existing Docker configurations +and scripts that reference `process_lidar.py` directly. + +Prefer using: python -m lidar_pipeline """ -import os -import sys -import json -import subprocess -from concurrent.futures import ProcessPoolExecutor, as_completed -from pathlib import Path -import numpy as np -import rasterio -from rasterio.transform import from_bounds -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -from matplotlib import rcParams -from scipy import ndimage -from scipy.stats import binned_statistic_2d -from tqdm import tqdm - -try: - import laspy -except ImportError as e: - print(f"Erreur import: {e}") - sys.exit(1) - -try: - from rasterio.warp import transform as warp_transform - HAS_WARP = True -except ImportError: - HAS_WARP = False - -# GPU acceleration via CuPy -try: - import cupy as cp - import cupyx.scipy.ndimage as cp_ndimage - HAS_GPU = True - _xp = cp # Default array module (GPU) - _gpu_info = cp.cuda.runtime.getDeviceProperties(0) - _gpu_name = _gpu_info['name'].decode() if isinstance(_gpu_info['name'], bytes) else str(_gpu_info['name']) - print(f"✓ GPU détectée: {_gpu_name}") - print(f" Mémoire GPU: {_gpu_info['totalGlobalMem'] // (1024**3)} Go") -except (ImportError, Exception): - HAS_GPU = False - _xp = np # Fallback CPU - - -def to_gpu(arr): - """Send array to GPU if available.""" - if HAS_GPU: - return cp.asarray(arr.astype(np.float64)) - return arr.astype(np.float64) - - -def to_cpu(arr): - """Bring array back to CPU (numpy).""" - if HAS_GPU and isinstance(arr, cp.ndarray): - return cp.asnumpy(arr) - return arr - - -def xp_gaussian_filter(arr, sigma): - """Gaussian filter using GPU if available.""" - if HAS_GPU and isinstance(arr, cp.ndarray): - return cp_ndimage.gaussian_filter(arr, sigma) - return ndimage.gaussian_filter(arr, sigma) - - -def xp_uniform_filter(arr, size): - """Uniform filter using GPU if available.""" - if HAS_GPU and isinstance(arr, cp.ndarray): - return cp_ndimage.uniform_filter(arr, size) - return ndimage.uniform_filter(arr, size) - - -def xp_minimum_filter(arr, footprint=None, size=None): - """Minimum filter using GPU if available.""" - if HAS_GPU and isinstance(arr, cp.ndarray): - return cp_ndimage.minimum_filter(arr, footprint=footprint, size=size) - return ndimage.minimum_filter(arr, footprint=footprint, size=size) - -rcParams['figure.dpi'] = 150 -rcParams['savefig.dpi'] = 300 -rcParams['font.size'] = 10 - - -class LidarArchaeoPipeline: - """Pipeline Python pur pour analyse archéologique LiDAR""" - - def __init__(self, input_dir, output_dir, resolution=0.5, workers=1): - self.input_dir = Path(input_dir) - self.output_dir = Path(output_dir) - self.resolution = resolution - self.workers = workers - self.temp_dir = self.output_dir / "temp" - - if not self.input_dir.exists(): - raise ValueError(f"Répertoire introuvable: {self.input_dir}") - - self.output_dir.mkdir(parents=True, exist_ok=True) - self.temp_dir.mkdir(exist_ok=True) - - 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]: - d.mkdir(exist_ok=True) - - print(f"✓ Pipeline initialisé (Python pur)") - print(f" Entrée : {self.input_dir}") - print(f" Sortie : {self.output_dir}") - print(f" Résolution : {resolution}m/px") - print(f" Workers : {workers}") - - def find_laz_files(self): - """Trouve tous les fichiers LAZ/LAS""" - files = list(self.input_dir.glob("*.laz")) + list(self.input_dir.glob("*.las")) - print(f"✓ {len(files)} fichier(s) LiDAR trouvé(s)") - return sorted(files) - - def check_tools(self): - """Vérifie que les outils sont disponibles""" - for name, cmd in [('pdal', 'pdal --version'), ('gdal', 'gdalinfo --version')]: - try: - subprocess.run(cmd.split(), capture_output=True, check=True) - print(f" ✓ {name}") - except (subprocess.CalledProcessError, FileNotFoundError): - return False - return True - - # ============ STEP 1: Classification ============ - - def create_pipeline_json(self, input_laz, output_las): - """Create PDAL pipeline for ground classification""" - pipeline = { - "pipeline": [ - str(input_laz), - { - "type": "filters.smrf", - "ignore": "Classification[7:7]", - "slope": 1.0, - "window": 16.0, - "threshold": 0.5, - "scalar": 1.25 - }, - { - "type": "filters.range", - "limits": "Classification[2:2]" - }, - { - "type": "writers.las", - "filename": str(output_las), - "extra_dims": "all" - } - ] - } - return json.dumps(pipeline) - - def classify_ground(self, laz_file): - """Classify ground points using SMRF filter""" - output_las = self.temp_dir / f"{laz_file.stem}_ground.las" - - if output_las.exists(): - print(f" ✓ Classification déjà effectuée") - return output_las - - pipeline_json = self.create_pipeline_json(laz_file, output_las) - pipeline_file = self.temp_dir / "pipeline.json" - - with open(pipeline_file, 'w') as f: - f.write(pipeline_json) - - try: - subprocess.run( - ["pdal", "pipeline", str(pipeline_file)], - capture_output=True, - check=True - ) - print(f" ✓ Classification sol terminée") - return output_las - except subprocess.CalledProcessError as e: - error_msg = e.stderr.decode() - print(f" ✗ Erreur PDAL: {error_msg}") - - # If error is about ReturnNumber=0, try filtering those points - if "ReturnNumber" in error_msg and "NumberOfReturns" in error_msg: - print(f" → Tentative de filtrage des points ReturnNumber=0...") - - # Use filters.range to keep only valid points (ReturnNumber >= 1) - filtered_pipeline = [ - { - "type": "readers.las", - "filename": str(laz_file) - }, - { - "type": "filters.range", - "limits": "ReturnNumber[1:],NumberOfReturns[1:]" - }, - { - "type": "filters.smrf", - "scalar": 1.25 - }, - { - "type": "filters.range", - "limits": "Classification[2:2]" - }, - { - "type": "writers.las", - "filename": str(output_las), - "extra_dims": "all" - } - ] - - filtered_json = json.dumps(filtered_pipeline) - with open(pipeline_file, 'w') as f: - f.write(filtered_json) - - try: - subprocess.run( - ["pdal", "pipeline", str(pipeline_file)], - capture_output=True, - check=True - ) - print(f" ✓ Classification sol terminée (points filtrés)") - return output_las - except subprocess.CalledProcessError as e2: - print(f" ✗ Erreur même avec filtrage: {e2.stderr.decode()}") - return None - - return None - - def run_pdal_pipeline(self, pipeline_json, output_file): - """Exécute un pipeline PDAL""" - pipeline_file = self.temp_dir / "temp_pipeline.json" - with open(pipeline_file, 'w') as f: - f.write(pipeline_json) - - try: - subprocess.run( - ["pdal", "pipeline", str(pipeline_file)], - capture_output=True, - check=True - ) - print(f" ✓ Pipeline terminé") - return output_file - except subprocess.CalledProcessError as e: - print(f" ✗ Erreur PDAL: {e.stderr.decode()}") - return None - - # ============ STEP 2: DTM Generation ============ - - def create_dtm_fast(self, las_file, basename): - """Create DTM using fast binning method with gap filling""" - print(f" → Génération DTM...") - - try: - import laspy - from scipy.ndimage import gaussian_filter - - # Read LAZ/LAS file - las = laspy.read(str(las_file)) - - # Get bounds - min_x, max_x = float(las.header.min[0]), float(las.header.max[0]) - min_y, max_y = float(las.header.min[1]), float(las.header.max[1]) - - # Calculate grid dimensions - width = int(np.ceil((max_x - min_x) / self.resolution)) - height = int(np.ceil((max_y - min_y) / self.resolution)) - - print(f" Bounds: X[{min_x:.1f}, {max_x:.1f}] Y[{min_y:.1f}, {max_y:.1f}]") - print(f" Grid: {width}x{height} pixels ({len(las.points):,} points)") - - # Fast binning method - print(f" Rasterisation rapide...") - - stat = binned_statistic_2d( - las.x, las.y, las.z, - statistic='mean', - bins=[width, height], - range=[[min_x, max_x], [min_y, max_y]] - ) - - dtm = stat.statistic.T - - # Flip Y axis so north is at the top (Y decreases from top to bottom in image) - dtm = dtm[::-1, :] - - # Fill gaps using interpolation - print(f" Remplissage des zones vides...") - - # Create mask of empty cells - mask = np.isnan(dtm) - - if np.any(mask): - from scipy.ndimage import distance_transform_edt - - # Calculate distance to nearest valid point for each cell - dist_to_nearest = distance_transform_edt(mask) - max_fill_distance = max(20, int(10 / self.resolution)) # Max 20px or 10m - - # Use nearest-neighbor interpolation for speed (memory efficient) - # but apply strong Gaussian smoothing to remove artifacts - from scipy import interpolate - - y_coords, x_coords = np.where(~mask) - z_values = dtm[~mask] - - if len(y_coords) > 100: - # Nearest neighbor interpolation (fast, memory-efficient) - interp = interpolate.NearestNDInterpolator( - np.column_stack((y_coords, x_coords)), - z_values - ) - - # Only interpolate where needed (memory-efficient) - y_missing, x_missing = np.where(mask) - dtm[y_missing, x_missing] = interp(y_missing, x_missing) - - # Apply Gaussian smoothing to entire DTM to remove blocky artifacts - # Stronger smoothing in areas that were interpolated - dtm_smooth = gaussian_filter(dtm, sigma=2.0) - - # Only use smoothed values for interpolated areas, keep original for valid data - # Also only fill within max distance from valid data - fill_mask = mask & (dist_to_nearest <= max_fill_distance) - dtm[fill_mask] = dtm_smooth[fill_mask] - - # For cells beyond max distance, leave as NaN - far_mask = mask & (dist_to_nearest > max_fill_distance) - dtm[far_mask] = np.nan - - # Save as GeoTIFF - output_tif = self.dtm_dir / f"{basename}_dtm.tif" - transform = from_bounds(min_x, min_y, max_x, max_y, width, height) - - with rasterio.open( - output_tif, - 'w', - driver='GTiff', - height=height, - width=width, - count=1, - dtype='float32', - crs='EPSG:2154', - transform=transform, - compress='lzw' - ) as dst: - dst.write(dtm.astype('float32'), 1) - - print(f" ✓ DTM créé: {output_tif.name}") - return output_tif - - except Exception as e: - print(f" ✗ Erreur DTM: {e}") - import traceback - traceback.print_exc() - return None - - # ============ STEP 3: Visualizations (Python) ============ - - def generate_hillshade(self, dem_file, basename): - """Generate hillshade using numpy""" - print(f" → Hillshade multidirectionnel...") - - output = self.vis_dir / f"{basename}_hillshade_multi.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Calculate gradients - dy, dx = np.gradient(dem) - - # Multi-directional hillshade (NW, NE, SW, SE) - azimuts = [315, 45, 225, 135] - altitude = 30 - hillshades = [] - - for az in azimuts: - az_rad = np.radians(az) - alt_rad = np.radians(altitude) - - # Calculate slope and aspect - slope = np.arctan(np.sqrt(dx**2 + dy**2)) - aspect = np.arctan2(dy, dx) - - # Hillshade formula - hs = np.sin(alt_rad) * np.sin(slope) + \ - np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect) - hillshades.append(np.clip(hs, 0, 1)) - - # Combine all directions - combined = np.mean(hillshades, axis=0) - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=combined.shape[0], - width=combined.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(combined.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur hillshade: {e}") - return None - - def generate_hillshade_ne(self, dem_file, basename): - """Generate hillshade from NW (azimuth 315deg) like topographic maps""" - print(f" → Hillshade Nord-Est (carte topo)...") - - output = self.vis_dir / f"{basename}_hillshade_ne.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Single-direction hillshade from NW (315°) - standard for topo maps - # Since DTM Y-axis is flipped (north at top, Y increases downward in image), - # use +dy in aspect calculation (not -dy) - azimuth = 315 - altitude = 45 - - dy, dx = np.gradient(dem) - slope = np.arctan(np.sqrt(dx**2 + dy**2)) - # Use +dy because DTM is already flipped (north at top) - aspect = np.arctan2(dy, dx) - - az_rad = np.radians(azimuth) - alt_rad = np.radians(altitude) - - hs = np.sin(alt_rad) * np.sin(slope) + \ - np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect) - hs = np.clip(hs, 0, 1) - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=hs.shape[0], - width=hs.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(hs.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur hillshade NE: {e}") - return None - - def generate_slope(self, dem_file, basename): - """Generate slope map""" - print(f" → Pente (Slope)...") - - output = self.vis_dir / f"{basename}_slope.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Calculate slope - dy, dx = np.gradient(dem) - slope = np.arctan(np.sqrt(dx**2 + dy**2)) * 180 / np.pi - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=slope.shape[0], - width=slope.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(slope.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur slope: {e}") - return None - - def generate_aspect(self, dem_file, basename): - """Generate aspect (orientation of slopes) map""" - print(f" → Aspect (Orientation)...") - - output = self.vis_dir / f"{basename}_aspect.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Calculate aspect (direction of slope) - dy, dx = np.gradient(dem) - aspect = np.arctan2(dy, dx) * 180 / np.pi - aspect = np.mod(aspect, 360) # Convert to 0-360 - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=aspect.shape[0], - width=aspect.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(aspect.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur aspect: {e}") - return None - - def generate_curvature(self, dem_file, basename): - """Generate curvature (terrain concavity/convexity) map""" - print(f" → Courbure (Curvature)...") - - output = self.vis_dir / f"{basename}_curvature.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Calculate curvature using Laplacian - dz_dx = np.gradient(dem, axis=1) - dz_dy = np.gradient(dem, axis=0) - d2z_dx2 = np.gradient(dz_dx, axis=1) - d2z_dy2 = np.gradient(dz_dy, axis=0) - - # Mean curvature - curvature = (d2z_dx2 + d2z_dy2) / 2 - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=curvature.shape[0], - width=curvature.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(curvature.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur curvature: {e}") - return None - - def generate_solar(self, dem_file, basename): - """Generate solar irradiance simulation""" - print(f" → Éclairage Solaire (Solar Irradiance)...") - - output = self.vis_dir / f"{basename}_solar.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Calculate gradients - dy, dx = np.gradient(dem) - - # Calculate slope and aspect - slope = np.arctan(np.sqrt(dx**2 + dy**2)) - aspect = np.arctan2(dy, dx) - - # Solar irradiance (morning sun - azimuth 90, altitude 30) - az_rad = np.radians(90) - alt_rad = np.radians(30) - - # Solar radiation formula - solar = np.sin(alt_rad) * np.sin(slope) + \ - np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect) - - # Clip negative values (shadows) - solar = np.clip(solar, 0, 1) - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=solar.shape[0], - width=solar.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(solar.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur solar: {e}") - return None - - def generate_lrm(self, dem_file, basename): - """Local Relief Model - deviation from local mean""" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → Local Relief Model{gpu_tag}...") - - output = self.vis_dir / f"{basename}_lrm.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - dem = to_gpu(dem_np) - local_mean = xp_gaussian_filter(dem, sigma=15/self.resolution) - lrm = dem - local_mean - lrm_np = to_cpu(lrm).astype(np.float32) - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=lrm_np.shape[0], - width=lrm_np.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(lrm_np, 1) - - return output - except Exception as e: - print(f" ✗ Erreur LRM: {e}") - return None - - def generate_svf(self, dem_file, basename): - """Sky-View Factor - ray-tracing on 16 azimuths. - For each pixel, trace rays in N directions, find the max horizon - angle in each direction, then SVF = (1/N) * sum(cos²(horizon_angle)). - Valleys/crevices have low SVF (obstructed sky), ridges/peaks have high SVF. - Uses GPU (CuPy) if available for acceleration.""" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → Sky-View Factor (ray-tracing){gpu_tag}...") - - output = self.vis_dir / f"{basename}_svf.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - rows, cols = dem_np.shape - res = self.resolution - - # Move to GPU if available - dem = to_gpu(dem_np) - - # 16 azimuth directions - n_dirs = 16 - angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx = np.cos(angles) - dy = np.sin(angles) - - max_dist = int(50 / res) # 50m search radius - - # Pad DEM with NaN for boundary handling - xp = cp if HAS_GPU else np - padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - - svf = xp.zeros_like(dem) - - for d_idx in range(n_dirs): - ddx, ddy = dx[d_idx], dy[d_idx] - horizon = xp.zeros_like(dem) - - for step in range(1, max_dist + 1): - px = int(round(ddx * step)) - py = int(round(ddy * step)) - dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) - if dist_m < res * 0.5: - continue - - elev_diff = padded[max_dist + py:max_dist + py + rows, - max_dist + px:max_dist + px + cols] - dem - - angle = xp.arctan2(elev_diff, dist_m) - horizon = xp.where(xp.isnan(angle), horizon, - xp.maximum(horizon, xp.nan_to_num(angle, nan=0))) - - svf += xp.cos(xp.pi / 2 - horizon) ** 2 - - svf /= n_dirs - - # Bring back to CPU for saving - svf_np = to_cpu(svf).astype(np.float32) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=svf_np.shape[0], - width=svf_np.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(svf_np, 1) - - return output - except Exception as e: - print(f" ✗ Erreur SVF: {e}") - return None - - return output - except Exception as e: - print(f" ✗ Erreur SVF: {e}") - return None - - def generate_openness(self, dem_file, basename, positive=True): - """Positive/Negative Openness - true zenith/nadir angle computation. - For each pixel, in 8 directions (N, NE, E, SE, S, SW, W, NW): - - Positive openness: max zenith angle (angle from vertical to highest visible terrain) - - Negative openness: max nadir angle (angle from vertical down to lowest terrain) - Result is averaged across all 8 directions. - Uses GPU (CuPy) if available for acceleration.""" - name = "positive_openness" if positive else "negative_openness" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → {name.replace('_', ' ').title()} (ray-tracing){gpu_tag}...") - - output = self.vis_dir / f"{basename}_{name}.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - rows, cols = dem_np.shape - res = self.resolution - - dem = to_gpu(dem_np) - xp = cp if HAS_GPU else np - - n_dirs = 8 - angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx = np.cos(angles) - dy = np.sin(angles) - - max_dist = int(50 / res) - - padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - - openness_sum = xp.zeros_like(dem) - - for d_idx in range(n_dirs): - ddx, ddy = dx[d_idx], dy[d_idx] - max_angle = xp.zeros_like(dem) - - for step in range(1, max_dist + 1): - px = int(round(ddx * step)) - py = int(round(ddy * step)) - dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) - if dist_m < res * 0.5: - continue - - elev_diff = padded[max_dist + py:max_dist + py + rows, - max_dist + px:max_dist + px + cols] - dem - - if positive: - angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) - else: - angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) - - max_angle = xp.where(xp.isnan(angle), max_angle, - xp.maximum(max_angle, xp.nan_to_num(angle, nan=0))) - - openness_sum += max_angle - - # Average across all directions (convert to degrees) - openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32) - with rasterio.open( - output, - 'w', - driver='GTiff', - height=openness_result.shape[0], - width=openness_result.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(openness_result.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur openness: {e}") - return None - - def generate_mslrm(self, dem_file, basename): - """Multi-Scale Relief Model (MSRM) - LRM computed at multiple scales - (sigma=5,10,25,50,100) and combined. Reveals features at all scales: - small (walls, ditches), medium (enclosures, tumulus), large (landscapes). - Uses GPU if available for acceleration.""" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...") - - output = self.vis_dir / f"{basename}_mslrm.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - dem = to_gpu(dem_np) - xp = cp if HAS_GPU else np - - # Compute LRM at multiple scales - sigmas = [5, 10, 25, 50, 100] - lrm_stack = [] - - for sigma in sigmas: - sigma_px = sigma / self.resolution - local_mean = xp_gaussian_filter(dem, sigma=sigma_px) - lrm = dem - local_mean - lrm_norm = lrm / max(xp.nanstd(lrm), 0.01) - lrm_stack.append(lrm_norm) - - # Combine: RMS of normalized LRM at all scales - mslrm = xp.sqrt(xp.mean(xp.array(lrm_stack) ** 2, axis=0)) - mslrm_np = to_cpu(mslrm).astype(np.float32) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=mslrm_np.shape[0], - width=mslrm_np.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(mslrm_np, 1) - - return output - except Exception as e: - print(f" ✗ Erreur MSRM: {e}") - return None - - def generate_tpi(self, dem_file, basename): - """Multi-Scale Topographic Position Index. - TPI = elevation - mean(neighborhood). - Computed at fine (5m) and broad (100m) scales to identify - ridges, valleys, and intermediate landforms. - Uses GPU if available for acceleration.""" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → TPI multi-echelle{gpu_tag}...") - - output = self.vis_dir / f"{basename}_tpi.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - dem = to_gpu(dem_np) - - fine_size = int(5 / self.resolution) - if fine_size % 2 == 0: - fine_size += 1 - tpi_fine = dem - xp_uniform_filter(dem, size=fine_size) - - broad_size = int(100 / self.resolution) - if broad_size % 2 == 0: - broad_size += 1 - tpi_broad = dem - xp_uniform_filter(dem, size=broad_size) - - # Combine: fine-scale weighted more for archaeological features - # Normalize each scale then weight: 0.6 fine + 0.4 broad - xp = cp if HAS_GPU else np - fine_std = max(float(xp.nanstd(tpi_fine)), 0.01) - broad_std = max(float(xp.nanstd(tpi_broad)), 0.01) - tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std) - tpi_np = to_cpu(tpi_combined).astype(np.float32) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=tpi_np.shape[0], - width=tpi_np.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(tpi_np, 1) - - return output - except Exception as e: - print(f" ✗ Erreur TPI: {e}") - return None - - def generate_vat(self, dem_file, basename): - """VAT composite (Visualisation for Archaeological Topography). - Fuses hillshade (R), slope (G), and SVF (B) into an RGB composite - that reveals all micro-topographic features in a single image.""" - print(f" → VAT composite...") - - output = self.vis_dir / f"{basename}_vat.tif" - - try: - # Generate required rasters if not already available - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - from scipy.ndimage import gaussian_filter - rows, cols = dem.shape - - # Hillshade (luminance) for R channel - gy, gx = np.gradient(dem, self.resolution) - slope = np.arctan(np.sqrt(gx**2 + gy**2)) - # Illumination from NW (azimuth 315°, altitude 45°) - alt_rad = np.radians(45) - az_rad = np.radians(315) - shading = (np.sin(alt_rad) * np.cos(slope) + - np.cos(alt_rad) * np.sin(slope) * - np.cos(az_rad - np.arctan2(gy, gx))) - hillshade = np.clip(shading, 0, 1) - - # Slope for G channel - slope_deg = np.degrees(slope) - slope_norm = np.clip(slope_deg / 45, 0, 1) - - # SVF (simplified fast version for composite) for B channel - n_dirs_svf = 8 - angles_svf = np.linspace(0, 2 * np.pi, n_dirs_svf, endpoint=False) - dx_svf = np.cos(angles_svf) - dy_svf = np.sin(angles_svf) - max_dist_svf = int(30 / self.resolution) - padded = np.pad(dem.astype(np.float64), max_dist_svf, mode='constant', - constant_values=np.nan) - svf = np.zeros_like(dem, dtype=np.float64) - for d_idx in range(n_dirs_svf): - horizon = np.zeros_like(dem, dtype=np.float64) - for step in range(1, max_dist_svf + 1): - px = int(round(dx_svf[d_idx] * step)) - py = int(round(dy_svf[d_idx] * step)) - dist_m = np.sqrt((dx_svf[d_idx] * step * self.resolution) ** 2 + - (dy_svf[d_idx] * step * self.resolution) ** 2) - if dist_m < self.resolution * 0.5: - continue - elev_diff = padded[max_dist_svf + py:max_dist_svf + py + rows, - max_dist_svf + px:max_dist_svf + px + cols] - dem - angle = np.arctan2(elev_diff, dist_m) - horizon = np.where(np.isnan(angle), horizon, - np.maximum(horizon, np.nan_to_num(angle, nan=0))) - svf += np.cos(np.pi / 2 - horizon) ** 2 - svf /= n_dirs_svf - - # Normalize each channel to 0-1 using percentile stretch - def stretch(arr, p_low=2, p_high=98): - valid = arr[~np.isnan(arr)] - if len(valid) == 0: - return arr - lo, hi = np.percentile(valid, (p_low, p_high)) - if hi - lo < 1e-10: - return np.zeros_like(arr) - return np.clip((arr - lo) / (hi - lo), 0, 1) - - r = stretch(hillshade, 1, 99) - g = stretch(slope_norm, 2, 95) - b = stretch(svf, 5, 95) - - # Stack as RGB (3-band GeoTIFF) - rgb = np.stack([r, g, b], axis=0).astype('float32') - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=rows, - width=cols, - count=3, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - for i in range(3): - dst.write(rgb[i], i + 1) - - return output - except Exception as e: - print(f" ✗ Erreur VAT: {e}") - return None - - def generate_depressions(self, dem_file, basename): - """Depression detection using hydrological sink filling. - Fills depressions in the DEM and computes the difference to identify - dolines, sinkholes, cavities, and enclosed basins.""" - print(f" → Detection depressions (hydrologique)...") - - output = self.vis_dir / f"{basename}_depressions.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Fill depressions using iterative approach (scipy) - # Priority-flood algorithm approximation - from scipy.ndimage import binary_dilation, generate_binary_structure - - # Replace NaN with a very high value for filling - dem_filled = dem.copy() - nodata_mask = np.isnan(dem_filled) - dem_filled[nodata_mask] = np.nanmax(dem) + 1000 - - # Iterative depression filling - # A pixel is a sink if it's lower than all its neighbors - # Fill by raising sinks to the minimum of their neighbors - changed = True - iterations = 0 - max_iter = 100 - - struct = generate_binary_structure(2, 2) # 8-connectivity - - while changed and iterations < max_iter: - # Find local minima (sinks) - from scipy.ndimage import minimum_filter - neighbor_min = minimum_filter(dem_filled, footprint=struct) - # A pixel is a sink if it's lower than the minimum of its neighbors - sinks = (dem_filled < neighbor_min) & ~nodata_mask - - if not np.any(sinks): - break - - # Fill sinks to the level of the lowest neighbor - new_dem = np.maximum(dem_filled, neighbor_min) - new_dem[nodata_mask] = np.nan - - changed = np.any(new_dem != dem_filled) - dem_filled = new_dem - iterations += 1 - - # Difference = depth of depressions - # Only keep positive differences (actual depressions) - depressions = dem_filled - dem - depressions[nodata_mask] = np.nan - # Only keep real depressions (positive depth) - depressions = np.where(depressions > 0.01, depressions, 0) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=depressions.shape[0], - width=depressions.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(depressions.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur depressions: {e}") - return None - - def generate_sailore(self, dem_file, basename): - """SAILORE - Self-Adaptive Improved Local Relief Model. - Kernel size adapts to local slope: flat areas get larger kernels, - steep areas get smaller kernels. Better than fixed LRM for heterogeneous terrain. - Uses GPU if available for acceleration.""" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → SAILORE (LRM adaptatif){gpu_tag}...") - - output = self.vis_dir / f"{basename}_sailore.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - dem = to_gpu(dem_np) - xp = cp if HAS_GPU else np - - gy, gx = xp.gradient(dem, self.resolution) - slope = xp.arctan(xp.sqrt(gx**2 + gy**2)) - slope_deg = xp.degrees(slope) - - sigma_min = 2.0 / self.resolution - sigma_max = 25.0 / self.resolution - slope_norm = xp.clip(slope_deg / 30.0, 0, 1) - adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min) - - lrm_fine = dem - xp_gaussian_filter(dem, sigma=sigma_min) - lrm_medium = dem - xp_gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2) - lrm_coarse = dem - xp_gaussian_filter(dem, sigma=sigma_max) - - w_fine = slope_norm - w_medium = 1 - 2 * xp.abs(slope_norm - 0.5) - w_coarse = 1 - slope_norm - w_total = w_fine + w_medium + w_coarse - w_total[w_total == 0] = 1 - - sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total - sailore_np = to_cpu(sailore).astype(np.float32) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=sailore_np.shape[0], - width=sailore_np.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(sailore_np, 1) - - return output - except Exception as e: - print(f" ✗ Erreur SAILORE: {e}") - return None - - def generate_geomorphons(self, dem_file, basename): - """Geomorphons - classification of terrain into 10 landform types - using ternary pattern of zenith/nadir angles in 8 directions. - Types: flat, peak, ridge, shoulder, spur, slope, hollow, footslope, - valley, pit.""" - print(f" → Geomorphons (10 formes de terrain)...") - - output = self.vis_dir / f"{basename}_geomorphons.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - rows, cols = dem.shape - - # Flatness threshold in meters - flat_threshold = 1.0 # 1m difference considered flat - - # 8 directions (N, NE, E, SE, S, SW, W, NW) - n_dirs = 8 - angles_g = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx_g = np.cos(angles_g) - dy_g = np.sin(angles_g) - - # Lookup distance in pixels - lookup = int(10 / self.resolution) # 10m lookup distance - - # Pad DEM - padded = np.pad(dem.astype(np.float64), lookup, mode='edge') - - # For each direction, determine if higher (+1), flat (0), or lower (-1) - pattern = np.zeros((n_dirs, rows, cols), dtype=np.int8) - - for d in range(n_dirs): - px = int(round(dx_g[d] * lookup)) - py = int(round(dy_g[d] * lookup)) - - neighbor = padded[lookup + py:lookup + py + rows, - lookup + px:lookup + px + cols] - - diff = neighbor - dem - pattern[d] = np.where(diff > flat_threshold, 1, - np.where(diff < -flat_threshold, -1, 0)) - - # Map ternary patterns to 10 landform types - # Using the standard geomorphon classification - # Count positive, negative, and flat directions - n_positive = np.sum(pattern == 1, axis=0) - n_negative = np.sum(pattern == -1, axis=0) - n_flat = np.sum(pattern == 0, axis=0) - - # Classification rules (standard geomorphons) - geomorphons = np.full((rows, cols), 0, dtype=np.int8) - # 0=flat, 1=peak, 2=ridge, 3=shoulder, 4=spur, - # 5=slope, 6=hollow, 7=footslope, 8=valley, 9=pit - - # Flat: all 8 directions are flat - geomorphons[n_flat == 8] = 0 - # Peak: all directions look up (all positive) - geomorphons[n_positive == 8] = 1 - # Pit: all directions look down (all negative) - geomorphons[n_negative == 8] = 9 - # Ridge: at least 3 positive, no negative, not all positive - geomorphons[(n_positive >= 3) & (n_negative == 0) & (n_positive < 8)] = 2 - # Valley: at least 3 negative, no positive, not all negative - geomorphons[(n_negative >= 3) & (n_positive == 0) & (n_negative < 8)] = 8 - # Shoulder: mixed positive+flat, some negative - geomorphons[(n_positive >= 2) & (n_negative >= 1) & (n_positive > n_negative)] = 3 - # Spur: some positive, more flat than negative - geomorphons[(n_positive >= 1) & (n_negative >= 1) & (n_positive > n_negative) & - (geomorphons == 0)] = 4 - # Slope: mostly flat with some negative or positive - geomorphons[(n_flat >= 4) & (n_positive + n_negative <= 4) & (geomorphons == 0)] = 5 - # Hollow: more negative than positive, some flat - geomorphons[(n_negative >= 1) & (n_positive >= 1) & (n_negative > n_positive) & - (geomorphons == 0)] = 6 - # Footslope: mixed negative+flat, some positive - geomorphons[(n_negative >= 2) & (n_positive >= 1) & (n_negative > n_positive) & - (geomorphons == 0)] = 7 - # Catch remaining as slope - geomorphons[geomorphons == 0] = 5 - - # Mask nodata - nodata_mask = np.isnan(dem) - geomorphons[nodata_mask] = 0 - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=rows, - width=cols, - count=1, - dtype='int8', - crs=crs, - transform=transform, - compress='lzw', - nodata=0 - ) as dst: - dst.write(geomorphons.astype('int8'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur Geomorphons: {e}") - return None - - def generate_roughness(self, dem_file, basename): - """Surface roughness - standard deviation of elevation in a window - plus roughness index (ratio of surface area to planar area). - Reveals anthropogenic surfaces (walls, paths) vs natural terrain.""" - print(f" → Rugosite de surface...") - - output = self.vis_dir / f"{basename}_roughness.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - from scipy.ndimage import generic_filter - - # Standard deviation in 5m window - window_size = int(5 / self.resolution) - if window_size % 2 == 0: - window_size += 1 - - # Compute std deviation (surface roughness) - def std_filter(arr): - return np.nanstd(arr) - - roughness = generic_filter(dem.astype(np.float64), std_filter, - size=window_size, mode='nearest') - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=roughness.shape[0], - width=roughness.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(roughness.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur rugosite: {e}") - return None - - def generate_anomalies(self, dem_file, basename): - """Statistical anomaly detection - z-score of local relief - plus Local Moran's I for spatial autocorrelation. - Identifies statistically significant topographic anomalies.""" - print(f" → Detection anomalies statistiques...") - - output = self.vis_dir / f"{basename}_anomalies.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Z-score of LRM (local relief) - from scipy.ndimage import gaussian_filter - - lrm = dem - gaussian_filter(dem, sigma=15 / self.resolution) - lrm_mean = np.nanmean(lrm) - lrm_std = max(np.nanstd(lrm), 0.01) - z_score = (lrm - lrm_mean) / lrm_std - - # Local Moran's I for spatial clustering - # I_i = (x_i - x_mean) * sum(w_ij * (x_j - x_mean)) / (S0 * sigma²) - # Simplified: use local mean of z-score as spatial lag - from scipy.ndimage import uniform_filter - - window = int(10 / self.resolution) - if window % 2 == 0: - window += 1 - - local_mean = uniform_filter(z_score, size=window) - local_var = uniform_filter(z_score ** 2, size=window) - local_mean ** 2 - local_var = np.maximum(local_var, 0.01) - - # Moran's I approximation: high positive = cluster, high negative = outlier - morans_i = z_score * (local_mean - np.nanmean(z_score)) / np.nanstd(z_score) - - # Combine: absolute z-score (anomaly magnitude) * sign of Moran's I (cluster type) - # High positive = significant elevated cluster, high negative = significant depression cluster - anomaly_score = np.abs(z_score) * np.sign(morans_i) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=anomaly_score.shape[0], - width=anomaly_score.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(anomaly_score.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur anomalies: {e}") - return None - - def generate_wavelet(self, dem_file, basename): - """Mexican Hat wavelet (Ricker wavelet) multi-scale analysis. - Continuous Wavelet Transform at multiple scales to detect - circular/circular features like tumulus, ring ditches, enclosures. - Uses GPU if available for acceleration.""" - gpu_tag = " [GPU]" if HAS_GPU else "" - print(f" → Ondelette Mexican Hat multi-echelle{gpu_tag}...") - - output = self.vis_dir / f"{basename}_wavelet.tif" - - try: - with rasterio.open(dem_file) as src: - dem_np = src.read(1) - transform = src.transform - crs = src.crs - - dem = to_gpu(dem_np) - xp = cp if HAS_GPU else np - - scales = [2, 5, 10, 20, 50] # meters - wavelet_stack = [] - - for scale_m in scales: - sigma_px = scale_m / self.resolution - # 2D Mexican Hat = Laplacian of Gaussian - if HAS_GPU: - from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace - response = -gpu_gaussian_laplace(dem, sigma=sigma_px) - else: - from scipy.ndimage import gaussian_laplace - response = -gaussian_laplace(to_cpu(dem).astype(np.float64), sigma=sigma_px) - response = to_gpu(response) - response /= max(float(xp.nanstd(response)), 0.01) - wavelet_stack.append(response) - - combined = xp.sqrt(xp.mean(xp.array(wavelet_stack) ** 2, axis=0)) - combined_np = to_cpu(combined).astype(np.float32) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=combined.shape[0], - width=combined.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(combined_np, 1) - - return output - except Exception as e: - print(f" ✗ Erreur ondelette: {e}") - return None - - def generate_texture(self, dem_file, basename): - """GLCM texture analysis on hillshade - contrast, entropy, homogeneity. - Detects anthropogenic surfaces (plowing, paths, walls) vs natural terrain.""" - print(f" → Texture GLCM...") - - output = self.vis_dir / f"{basename}_texture.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # First create a hillshade for texture analysis - gy, gx = np.gradient(dem, self.resolution) - slope = np.arctan(np.sqrt(gx**2 + gy**2)) - alt_rad = np.radians(45) - az_rad = np.radians(315) - shading = (np.sin(alt_rad) * np.cos(slope) + - np.cos(alt_rad) * np.sin(slope) * - np.cos(az_rad - np.arctan2(gy, gx))) - hillshade = np.clip(shading, 0, 1) - - # Quantize to 256 levels for GLCM - valid = hillshade[~np.isnan(hillshade)] - if len(valid) == 0: - raise ValueError("No valid data for texture analysis") - lo, hi = np.percentile(valid, (1, 99)) - img = np.clip((hillshade - lo) / max(hi - lo, 0.001), 0, 1) - img_uint8 = (img * 255).astype(np.uint8) - - # Compute GLCM texture using sklearn-style approach - # Simplified GLCM: local variance (proxy for contrast) and local entropy - from scipy.ndimage import generic_filter - - window = int(5 / self.resolution) - if window % 2 == 0: - window += 1 - - # Local variance (contrast proxy) - def local_variance(arr): - arr_f = arr.astype(np.float64) - return np.var(arr_f) - - contrast = generic_filter(img_uint8.astype(np.float64), local_variance, - size=window, mode='nearest') - - # Local entropy approximation - def local_entropy(arr): - arr_f = arr.astype(np.float64) - hist, _ = np.histogram(arr_f, bins=16, range=(0, 256)) - hist = hist / max(hist.sum(), 1) - hist = hist[hist > 0] - return -np.sum(hist * np.log2(hist)) - - entropy = generic_filter(img_uint8.astype(np.float64), local_entropy, - size=window, mode='nearest') - - # Local homogeneity (inverse difference moment proxy) - def local_homogeneity(arr): - arr_f = arr.astype(np.float64) - mean_val = np.mean(arr_f) - return np.mean(1.0 / (1.0 + (arr_f - mean_val) ** 2)) - - homogeneity = generic_filter(img_uint8.astype(np.float64), local_homogeneity, - size=window, mode='nearest') - - # Combine into composite texture index - # Normalize each component - def norm(arr): - valid = arr[~np.isnan(arr)] - if len(valid) == 0: - return arr - std = max(np.std(valid), 0.01) - return (arr - np.mean(valid)) / std - - texture_combined = (0.4 * norm(contrast) + - 0.4 * norm(entropy) - - 0.2 * norm(homogeneity)) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=texture_combined.shape[0], - width=texture_combined.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(texture_combined.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur texture GLCM: {e}") - return None - - def generate_flow(self, dem_file, basename): - """Flow accumulation using D8 algorithm. - Identifies drainage patterns, ditches, and enclosure boundaries. - Uses scipy for efficient pit-filling and flow routing.""" - print(f" → Accumulation de flux D8...") - - output = self.vis_dir / f"{basename}_flow.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - rows, cols = dem.shape - nodata_mask = np.isnan(dem) - - # Fill pits first for proper flow routing - from scipy.ndimage import minimum_filter, generate_binary_structure - - # Simple pit-filling: iterative depression filling - dem_filled = dem.copy() - dem_filled[nodata_mask] = np.nanmax(dem) + 1000 - - struct = generate_binary_structure(2, 2) # 8-connectivity - for _ in range(50): - neighbor_min = minimum_filter(dem_filled, footprint=struct) - sinks = (dem_filled < neighbor_min) & ~nodata_mask - if not np.any(sinks): - break - dem_filled = np.where(sinks, neighbor_min, dem_filled) - - dem_filled[nodata_mask] = np.nan - - # D8 flow direction: each cell flows to steepest descent neighbor - # 8 neighbors: E, SE, S, SW, W, NW, N, NE - dx8 = [1, 1, 0, -1, -1, -1, 0, 1] - dy8 = [0, 1, 1, 1, 0, -1, -1, -1] - # Distances (cardinal=1, diagonal=sqrt(2)) - dist8 = [1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2), 1.0, np.sqrt(2)] - - # Compute flow direction for each cell - flow_dir = np.full((rows, cols), -1, dtype=np.int8) - max_slope = np.full((rows, cols), 0.0) - - padded = np.pad(dem_filled, 1, mode='constant', constant_values=np.nanmax(dem_filled) + 10000) - - for d in range(8): - # Neighbor elevation - nx = 1 + dx8[d] - ny = 1 + dy8[d] - neighbor_elev = padded[ny:ny + rows, nx:nx + cols] - # Slope to neighbor (positive = downhill) - slope = (dem_filled - neighbor_elev) / (dist8[d] * self.resolution) - slope[nodata_mask] = -1 - - # Update flow direction if this neighbor has steepest descent - better = slope > max_slope - flow_dir[better] = d - max_slope[better] = slope[better] - - # Flow accumulation: count upstream cells - # Process cells from highest to lowest elevation - flat_dem = dem_filled[~nodata_mask].flatten() - valid_indices = np.where(~nodata_mask.flatten())[0] - # Sort by elevation (highest first = process upstream first) - sort_order = valid_indices[np.argsort(-flat_dem)] - - flow_acc = np.ones((rows, cols), dtype=np.float32) - flow_acc[nodata_mask] = 0 - - # Accumulate flow - for idx in sort_order: - r, c = divmod(idx, cols) - d = flow_dir[r, c] - if d < 0: - continue - # Downstream cell - nr, nc = r + dy8[d], c + dx8[d] - if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]: - flow_acc[nr, nc] += flow_acc[r, c] - - # Log-scale for visualization - flow_log = np.log1p(flow_acc) - - with rasterio.open( - output, - 'w', - driver='GTiff', - height=flow_log.shape[0], - width=flow_log.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(flow_log.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur flux: {e}") - return None - """Detect cavities/depressions below ground level. - Uses local relief (difference between smoothed ground and actual terrain) - to highlight depressions, dolines, sinkholes and any below-ground features. - Surface features (mounds, walls) are suppressed so only depressions show.""" - print(f" → Detection cavites (gouffres)...") - - output = self.vis_dir / f"{basename}_cavity.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - from scipy.ndimage import uniform_filter, minimum_filter, gaussian_filter - from scipy.stats import binned_statistic_2d - - # Method 1: Local Relief - compare terrain to smoothed reference - # Large window smoothing captures the regional ground level - # Small depressions (cavities, dolines) appear as negative values - window_large = int(100 / self.resolution) # 100m window - if window_large % 2 == 0: - window_large += 1 - - # Regional ground level (smoothed over large area) - regional_ground = uniform_filter(dem, size=window_large) - - # Local relief = actual terrain minus regional ground - # Positive = above regional ground (hills, walls) - # Negative = below regional ground (cavities, dolines, sinkholes) - local_relief = dem - regional_ground - - # Method 2: Minimum filter to detect sharp local depressions - # Compares each point to its minimum neighborhood - window_small = int(20 / self.resolution) | 1 # 20m window - local_min = minimum_filter(dem, size=window_small) - - # Difference between local minimum and the actual terrain - # If terrain is much higher than local minimum, there's a deep pit nearby - pit_depth = dem - local_min - # Only keep negative values (terrain lower than neighbors = depression) - # This detects sharp, small depressions - - # Combine both methods: - # Method 1 catches large gentle depressions (dolines) - # Method 2 catches sharp small depressions (sinkholes, wells) - # Take the most negative (deepest) signal - combined = np.minimum(local_relief, -pit_depth) - - # Only show depressions (negative values = below ground) - # Set positive values (hills, walls) to 0 so they don't obscure cavities - combined = np.where(combined < 0, -combined, 0) - - # Slight smoothing to reduce noise - combined = gaussian_filter(combined, sigma=1.0) - - # Mask NaN areas - combined = np.where(np.isnan(dem), np.nan, combined) - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=combined.shape[0], - width=combined.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(combined.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur cavity detection: {e}") - import traceback - traceback.print_exc() - return None - - def generate_nodata_mask(self, dem_file, basename): - """Generate a map showing areas without LiDAR coverage""" - print(f" → Carte couverture LiDAR...") - - output = self.vis_dir / f"{basename}_nodata.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - - # Create binary mask: 1 = has data, 0 = no data (NaN) - # Also create distance from data edge to show interpolated areas - from scipy.ndimage import distance_transform_edt - - has_data = (~np.isnan(dem)).astype(np.float32) - no_data = np.isnan(dem).astype(np.float32) - - # Distance to nearest data point (in pixels) - if np.any(no_data.astype(bool)): - dist_from_data = distance_transform_edt(no_data.astype(bool)) - else: - dist_from_data = np.zeros_like(dem) - - # Normalize: 0 = has data, increasing values = further from data - # Cap at a reasonable distance - max_dist = 50 # pixels - coverage = np.clip(dist_from_data / max_dist, 0, 1) - - # Areas with data = 0, areas without data > 0 - # For visualization, use negative values for data and positive for gaps - result = np.where(np.isnan(dem), coverage, -0.1) # -0.1 for data areas - - # Save - with rasterio.open( - output, - 'w', - driver='GTiff', - height=result.shape[0], - width=result.shape[1], - count=1, - dtype='float32', - crs=crs, - transform=transform, - compress='lzw' - ) as dst: - dst.write(result.astype('float32'), 1) - - return output - except Exception as e: - print(f" ✗ Erreur nodata mask: {e}") - import traceback - traceback.print_exc() - return None - - def download_ign_tiles(self, min_x, max_x, min_y, max_y, layer, zoom_level=15): - """Download IGN WMTS tiles for the given bounds using Web Mercator (PM). - Converts Lambert 93 bounds to lat/lon then to Web Mercator tile coordinates.""" - import urllib.request - import io - import math - from PIL import Image as PILImage - - # Convert Lambert 93 bounds to WGS84 lat/lon - try: - l93_xs = [min_x, max_x] - l93_ys = [max_y, min_y] # NW, SE - lons, lats = warp_transform('EPSG:2154', 'EPSG:4326', l93_xs, l93_ys) - nw_lat, nw_lon = lats[0], lons[0] - se_lat, se_lon = lats[1], lons[1] - except Exception: - print(f" ✗ Impossible de convertir les coordonnees") - return None - - # IGN WMTS endpoint (free, no API key needed) - wmts_url = "https://data.geopf.fr/wmts" - tile_matrix_set = "PM" # Web Mercator (native cache) - tile_size = 256 - - # Web Mercator resolution at each zoom level (meters per pixel at equator) - # and tile calculation - n = 2 ** zoom_level # Number of tiles at this zoom level - - # Calculate tile range from lat/lon bounds - # Web Mercator formulas - def lat_lon_to_tile(lat, lon, zoom): - n = 2 ** zoom - col = int((lon + 180) / 360 * n) - lat_rad = math.radians(lat) - row = int((1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n) - return col, row - - col_min, row_min = lat_lon_to_tile(nw_lat, nw_lon, zoom_level) - col_max, row_max = lat_lon_to_tile(se_lat, se_lon, zoom_level) - - # Calculate Web Mercator pixel coordinates for compositing - def lat_lon_to_px(lat, lon, zoom): - n = 2 ** zoom - px_x = (lon + 180) / 360 * n * tile_size - lat_rad = math.radians(lat) - px_y = (1 - math.log(math.tan(lat_rad) + 1 / math.cos(lat_rad)) / math.pi) / 2 * n * tile_size - return px_x, px_y - - nw_px_x, nw_px_y = lat_lon_to_px(nw_lat, nw_lon, zoom_level) - se_px_x, se_px_y = lat_lon_to_px(se_lat, se_lon, zoom_level) - - out_width = int(se_px_x - nw_px_x) - out_height = int(se_px_y - nw_px_y) - - if out_width <= 0 or out_height <= 0 or out_width > 5000 or out_height > 5000: - return None - - composite = np.full((out_height, out_width, 3), 255, dtype=np.uint8) - - # Download and assemble tiles - tiles_downloaded = 0 - fmt = "image/png" if 'PLAN' in layer else "image/jpeg" - - for col in range(col_min, col_max + 1): - for row in range(row_min, row_max + 1): - url = ( - f"{wmts_url}?SERVICE=WMTS&VERSION=1.0.0&REQUEST=GetTile" - f"&LAYER={layer}&STYLE=normal" - f"&TILEMATRIXSET={tile_matrix_set}" - f"&TILEMATRIX={zoom_level}&TILECOL={col}&TILEROW={row}" - f"&FORMAT={fmt}" - ) - - try: - req = urllib.request.Request(url, headers={'User-Agent': 'Mozilla/5.0'}) - with urllib.request.urlopen(req, timeout=10) as response: - tile_data = response.read() - tile_img = PILImage.open(io.BytesIO(tile_data)).convert('RGB') - tile_arr = np.array(tile_img) - - # Calculate tile position in Web Mercator pixels - tile_origin_x = col * tile_size - tile_origin_y = row * tile_size - - # Position in composite image (relative to NW corner) - px_x = int(tile_origin_x - nw_px_x) - px_y = int(tile_origin_y - nw_px_y) - - # Paste tile into composite - x_off = max(0, -px_x) - y_off = max(0, -px_y) - dst_x_start = max(0, px_x) - dst_y_start = max(0, px_y) - dst_x_end = min(out_width, px_x + tile_size) - dst_y_end = min(out_height, px_y + tile_size) - - src_x = x_off - src_y = y_off - src_w = dst_x_end - dst_x_start - src_h = dst_y_end - dst_y_start - - if src_w > 0 and src_h > 0 and tile_arr.shape[0] >= src_y + src_h and tile_arr.shape[1] >= src_x + src_w: - composite[dst_y_start:dst_y_end, dst_x_start:dst_x_end] = \ - tile_arr[src_y:src_y+src_h, src_x:src_x+src_w] - tiles_downloaded += 1 - - except Exception as e: - if tiles_downloaded == 0 and col == col_min and row == row_min: - print(f" ✗ Erreur tuile IGN: {e}") - continue - - print(f" → {tiles_downloaded} tuiles IGN telechargees ({layer})") - if tiles_downloaded == 0: - return None - return composite - - def generate_ign_overlay(self, dem_file, basename, layer, title, legend_label, description, out_suffix): - """Generate an IGN basemap overlay visualization""" - print(f" → {title}...") - - output = self.vis_dir / f"{basename}_{out_suffix}.tif" - - try: - with rasterio.open(dem_file) as src: - dem = src.read(1) - transform = src.transform - crs = src.crs - height, width = dem.shape - - top_left_x = transform.c - top_left_y = transform.f - pixel_size_x = transform.a - pixel_size_y = abs(transform.e) - - min_x = top_left_x - max_x = top_left_x + width * pixel_size_x - max_y = top_left_y - min_y = top_left_y - height * pixel_size_y - - # Choose zoom level based on resolution - # At 0.5m/px, use zoom level 15 (3.19m/px IGN) - zoom = 15 - - result = self.download_ign_tiles(min_x, max_x, min_y, max_y, layer, zoom) - - if result is None: - print(f" ✗ Impossible de telecharger les tuiles IGN") - return None - - ign_img = result - - # Resize IGN image to match our DTM dimensions - from PIL import Image as PILImage - ign_pil = PILImage.fromarray(ign_img) - ign_resized = ign_pil.resize((width, height), PILImage.LANCZOS) - ign_arr = np.array(ign_resized) - - # Save as GeoTIFF - with rasterio.open( - output, - 'w', - driver='GTiff', - height=height, - width=width, - count=3, - dtype='uint8', - crs=crs, - transform=transform, - ) as dst: - for i in range(3): - dst.write(ign_arr[:, :, i], i + 1) - - return output - - except Exception as e: - print(f" ✗ Erreur IGN overlay: {e}") - import traceback - traceback.print_exc() - return None - - # ============ STEP 4: PNG Conversion ============ - - def tif_to_png(self, tif_file): - """Convert GeoTIFF to visualization PNG with GPS coordinates, external legend/scale""" - if not tif_file or not tif_file.exists(): - return None - - jpg_file = self.vis_dir / f"{tif_file.stem}.webp" - - try: - with rasterio.open(tif_file) as src: - # Check if this is a 3-band RGB image (ortho/topo) - is_rgb = src.count >= 3 and ('ortho' in str(tif_file) or 'topo' in str(tif_file)) - - if is_rgb: - # Read all 3 bands for RGB - data = src.read([1, 2, 3]) # Shape: (3, height, width) - data = np.moveaxis(data, 0, -1) # Shape: (height, width, 3) - else: - data = src.read(1) - - nodata = src.nodata - transform = src.transform - crs = src.crs - - # Get geographic bounds from transform - if is_rgb: - height, width, _ = data.shape - else: - height, width = data.shape - top_left_x = transform.c - top_left_y = transform.f - pixel_size_x = transform.a - pixel_size_y = abs(transform.e) - - min_x = top_left_x - max_x = top_left_x + width * pixel_size_x - max_y = top_left_y - min_y = top_left_y - height * pixel_size_y - - # Convert Lambert 93 corners to GPS (WGS84) - gps_coords = {} - if HAS_WARP and crs is not None: - try: - # Corner coordinates in Lambert 93 - l93_xs = [min_x, max_x, min_x, max_x] - l93_ys = [max_y, max_y, min_y, min_y] - lons, lats = warp_transform(crs, 'EPSG:4326', l93_xs, l93_ys) - gps_coords = { - 'NW': (lats[0], lons[0]), - 'NE': (lats[1], lons[1]), - 'SW': (lats[2], lons[2]), - 'SE': (lats[3], lons[3]), - } - # Also compute tick positions along edges - n_ticks = 5 - # Bottom edge (X axis): longitude ticks - tick_l93_x = np.linspace(min_x, max_x, n_ticks) - tick_l93_y_bottom = np.full(n_ticks, min_y) - tick_lons, tick_lats = warp_transform(crs, 'EPSG:4326', tick_l93_x, tick_l93_y_bottom) - gps_coords['x_ticks'] = list(zip(tick_lons, tick_lats)) - # Left edge (Y axis): latitude ticks - tick_l93_y = np.linspace(min_y, max_y, n_ticks) - tick_l93_x_left = np.full(n_ticks, min_x) - tick_lons_y, tick_lats_y = warp_transform(crs, 'EPSG:4326', tick_l93_x_left, tick_l93_y) - gps_coords['y_ticks'] = list(zip(tick_lons_y, tick_lats_y)) - except Exception: - gps_coords = {} - - if nodata is not None and not is_rgb: - data = np.ma.masked_where((data == nodata) | np.isnan(data), data) - - # For RGB images (ortho/topo), skip normalization and display directly - if is_rgb: - # RGB image: skip all normalization, display as-is - pass - else: - # Get valid data for normalization - valid_data = data.compressed() if hasattr(data, 'compressed') else data.flatten() - valid_data = valid_data[~np.isnan(valid_data)] - - # Determine colormap, normalization, and legend info - if is_rgb: - # RGB images: no colormap needed, display directly - if 'ortho' in str(tif_file): - title = "Photographie Aerienne IGN" - legend_label = "Orthophotographie\nImage aerienne" - description = "Photographie aerienne IGN" - elif 'topo' in str(tif_file): - title = "Carte Topographique IGN" - legend_label = "Carte IGN\nPlan topographique" - description = "Carte topographique IGN (Plan IGN)" - else: - title = tif_file.stem.replace('_', ' ').title() - legend_label = "Image RGB" - description = "" - elif 'hillshade' in str(tif_file): - cmap = 'gray' - vmin, vmax = 0, 1 - data = np.clip(data, vmin, vmax) - title = "Hillshade Multidirectionnel" - legend_label = "Illumination combinee de 6 directions\nBlanc = Face eclairee | Noir = Zone d'ombre" - description = "Ombres portees revelant micro-relief (murs, fosses, terrasses)" - elif 'slope' in str(tif_file): - cmap = 'inferno' - vmin, vmax = 0, np.percentile(valid_data, 95) - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Pente (Inclinaison du terrain)" - legend_label = f"Inclinaison en degres\nMin: {vmin:.1f} | Max: {vmax:.1f}\nClair = Forte pente | Sombre = Terrain plat" - description = "Murs, talus et bords ressortent en clair - terrain plat en sombre" - elif 'aspect' in str(tif_file): - cmap = 'hsv' - vmin, vmax = 0, 360 - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Aspect (Direction des pentes)" - legend_label = "Direction vers laquelle le terrain descend\nRouge=Nord | Vert=Est | Cyan=Sud | Bleu=Ouest" - description = "Orientation des pentes - utile pour distinguer structures des formes naturelles" - elif 'curvature' in str(tif_file): - cmap = 'RdYlBu_r' - vmax = max(abs(np.percentile(valid_data, 5)), abs(np.percentile(valid_data, 95)), 0.001) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Courbure (Convexite/Concavite du terrain)" - legend_label = f"Changement de pente (1/m)\nRouge = Convexe (sommet de mur, levee)\nBleu = Concave (fond de fosse, depression)" - description = "Detecte les ruptures de pente - utile pour bords de terrasses et levees" - elif 'svf' in str(tif_file): - cmap = 'viridis' - vmin, vmax = np.percentile(valid_data, (5, 95)) - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Sky-View Factor (Ray-tracing 16 directions)" - legend_label = "Fraction de ciel visible depuis chaque point\nSombre = Encaisse (fosses, vallees, rues)\nClair = Degage (sommets, plateformes, plateaux)" - description = "Ray-tracing sur 16 azimuts - elimine l'eclairage, detecte structures lineaires et enclos" - elif 'mslrm' in str(tif_file): - cmap = 'RdBu_r' - vmax = max(abs(np.percentile(valid_data, 2)), abs(np.percentile(valid_data, 98)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "MSRM - Multi-Scale Relief Model (5 echelles)" - legend_label = f"Relief combine a 5 echelles (5m a 100m)\nRouge = Surelevation (mur, tumulus, levee)\nBleu = Depression (fosse, douve, fosse)\n\nDifference avec LRM:\nLRM = 1 echelle (15m)\nMSRM = 5 echelles combinees\nMSRM detecte du micro au macro" - description = "Combine LRM a 5 echelles - detecte structures de 5m a 100m simultanement" - elif 'lrm' in str(tif_file): - cmap = 'RdBu_r' - vmax = max(abs(np.percentile(valid_data, 2)), abs(np.percentile(valid_data, 98)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "LRM - Local Relief Model (echelle unique 15m)" - legend_label = f"Ecart local par rapport au terrain moyen (m)\nRouge = Surelevation (+{vmax:.2f}m)\nBleu = Depression ({vmin:.2f}m)\nNoyau gaussien unique de 15m" - description = "Micro-relief a 15m seulement - voir MSRM pour toutes les echelles" - elif 'positive' in str(tif_file): - cmap = 'YlOrBr' - vmin, vmax = np.percentile(valid_data, (10, 98)) - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Openness Positive (ouverture vers le haut)" - legend_label = f"Angle d'ouverture vers le haut (deg)\nClair = Vue degagee vers le ciel (sommets, plateaux)\nSombre = Vue bouchee (vallees encaissees)" - description = "Ray-tracing 8 directions - complementaire de la negative pour detecter cretes" - elif 'negative' in str(tif_file): - cmap = 'GnBu_r' - vmin, vmax = np.percentile(valid_data, (10, 98)) - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Openness Negative (ouverture vers le bas)" - legend_label = f"Angle d'ouverture vers le bas (deg)\nClair = Surplomb (bords de fosse, grottes)\nSombre = Terrain plat (fonds de vallee)\nMeilleur detecteur de cavites et dolines" - description = "Ray-tracing 8 directions - detecte fosses, dolines, souterrains, dolines" - elif 'tpi' in str(tif_file): - cmap = 'BrBG' - vmax = max(abs(np.percentile(valid_data, 2)), abs(np.percentile(valid_data, 98)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "TPI - Topographic Position Index (2 echelles)" - legend_label = "Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fosse, vallee)\nVert/Clair = Plus haut que le voisinage (crete, plateau)\nCombine echelle fine 5m + large 100m" - description = "Identifie la position topographique - utile pour repeter cretes vs vallees a grande echelle" - elif 'depressions' in str(tif_file): - cmap = 'YlOrRd' - vmin, vmax = 0, max(np.percentile(valid_data, 99), 0.1) - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Depressions (Remplissage hydrologique)" - legend_label = f"Profondeur des cuvettes detectees (m)\nTransparent = Pas de depression\nJaune = Depression legere | Rouge = Depression profonde\nMax: {vmax:.2f}m" - description = "Simule le remplissage d'eau - detecte dolines, sinkholes, cuvettes et zones inondables" - elif 'sailore' in str(tif_file): - cmap = 'seismic' - vmax = max(abs(np.percentile(valid_data, 2)), abs(np.percentile(valid_data, 98)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "SAILORE - LRM Auto-Adaptatif" - legend_label = f"Relief local adaptatif (m)\nRouge = Surelevation | Bleu = Depression\n\nDifference avec LRM/MSRM:\nLRM = noyau fixe 15m\nMSRM = 5 noyaux fixes\nSAILORE = noyau adapte a la pente\nPlat=grand noyau | Pente=petit noyau" - description = "Noyau qui s'adapte a la pente locale - terrain plat=grand noyau, pente=petit noyau" - elif 'roughness' in str(tif_file): - cmap = 'magma' - vmin, vmax = 0, np.percentile(valid_data, 97) - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Rugosite de Surface (Ecart-type local 5m)" - legend_label = f"Irregularite du terrain dans un voisinage de 5m\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (vegetation, ruines, pierres)\nMax: {vmax:.2f}m" - description = "Mesure la variabilite locale - surfaces anthropiques lisses vs naturelles rugueuses" - elif 'anomalies' in str(tif_file): - cmap = 'coolwarm' - vmax = max(abs(np.percentile(valid_data, 5)), abs(np.percentile(valid_data, 95)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Anomalies Statistiques (Z-score x Moran's I)" - legend_label = "Anomalies topographiques significatives\nRouge vif = Surelevation anormale (mur, tumulus)\nBleu vif = Depression anormale (fosse, doline)\nBlanc/gris = Normal\n\nCombine z-score (intensite) et\nMoran's I (regroupement spatial)" - description = "Detecte uniquement les anomalies statistiquement significatives - filtre le bruit de fond" - elif 'wavelet' in str(tif_file): - cmap = 'cividis' - vmax = max(abs(np.percentile(valid_data, 2)), abs(np.percentile(valid_data, 98)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Ondelette Mexican Hat (CWT multi-echelle)" - legend_label = "Reponse de la transformee en ondelette a 5 echelles\nEchelles: 2m, 5m, 10m, 20m, 50m\n\nClair = Structure detectee a cette echelle\nSombre = Pas de structure\n\nOptimise pour formes circulaires:\ntumulus, enclos, fosses annulaires" - description = "Transformee en ondelette 2D - excellente pour detecter structures circulaires" - elif 'texture' in str(tif_file): - cmap = 'inferno' - vmax = max(abs(np.percentile(valid_data, 2)), abs(np.percentile(valid_data, 98)), 0.1) - vmin, vmax = -vmax, vmax - data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Texture GLCM (Contraste + Entropie - Homogeneite)" - legend_label = "Analyse de la texture du relief (fenetre 5m)\nClair = Texture heterogene (labour, ruines, sol perturbe)\nSombre = Texture homogene (sol nu, route, zone plate)\n\nCombine contraste, entropie et homogeneite" - description = "Distingue surfaces anthropiques (labour, chemins) des naturelles" - elif 'flow' in str(tif_file): - cmap = 'Blues' - vmin, vmax = 0, np.percentile(valid_data, 98) - data = np.clip((data - vmin) / max(vmax - vmin, 0.01), 0, 1) - title = "Accumulation de Flux Hydrologique (D8)" - legend_label = f"Logarithme de l'accumulation d'eau\nBlanc = Pas de collecte (sommet, crete)\nBleu fonce = Collecte maximale (thalweg, fosse)\n\nSimule l'ecoulement de l'eau de pluie\nDetecte fosses d'enceinte, routes et drainage" - description = "Algorithme D8 - simule le cheminement de l'eau pour detecter fosses et routes antiques" - elif 'ortho' in str(tif_file): - # Aerial photo - display RGB directly - title = "Photographie Aerienne IGN" - legend_label = "Orthophotographie\nImage aerienne" - description = "Photographie aerienne IGN" - # For 3-band RGB, display directly - if data.ndim == 3: - # Already RGB, just display - pass - else: - cmap = 'gray' - # Skip normalization for RGB - elif 'topo' in str(tif_file): - # Topographic map - display RGB directly - title = "Carte Topographique IGN" - legend_label = "Carte IGN\nPlan topographique" - description = "Carte topographique IGN (Plan IGN)" - # For 3-band RGB, display directly - if data.ndim == 3: - pass - else: - cmap = 'gray' - else: - cmap = 'terrain' - p2, p98 = np.percentile(valid_data, (2, 98)) - data = np.clip((data - p2) / (p98 - p2), 0, 1) - title = tif_file.stem.replace('_', ' ').title() - legend_label = "Altitude normalisee" - description = "" - - # ---- Create figure with external margins ---- - # Layout: map with colorbar, info bar below - # Map aspect ratio determines the figure height dynamically - from matplotlib.gridspec import GridSpec - # Keep image width-based sizing, add room for colorbar and bottom info - fig_width = 20 - # Map height based on data aspect ratio, plus margins for title and info - map_aspect = height / width - fig = plt.figure(figsize=(fig_width, fig_width * map_aspect * 0.7 + 2.5), - facecolor='white') - gs = GridSpec(2, 1, height_ratios=[1.0, 0.06], - hspace=0.04, figure=fig, - left=0.06, right=0.88, top=0.93, bottom=0.08) - - # Map row (with title above) - ax = fig.add_subplot(gs[0]) - # Display data - RGB or single-band - if is_rgb: - im = ax.imshow(data, aspect='equal', origin='upper') - else: - im = ax.imshow(data, cmap=cmap, aspect='equal', origin='upper') - - # Title above the map - ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10) - - # Colorbar on the right side (not for RGB images) - if not is_rgb: - cbar = plt.colorbar(im, ax=ax, pad=0.02, shrink=0.85, aspect=30) - cbar.ax.tick_params(labelsize=9, width=1.5) - cbar.outline.set_linewidth(1.5) - cbar.set_label(legend_label, fontsize=10, fontweight='bold') - else: - # For other RGB images (ortho/topo), add legend label as text on the right - ax.text(1.02, 0.5, legend_label, transform=ax.transAxes, - fontsize=10, fontweight='bold', rotation=90, - verticalalignment='center', horizontalalignment='left') - - # ---- GPS coordinate tick labels ---- - if gps_coords and 'x_ticks' in gps_coords: - # X-axis: longitude along bottom - x_pixel_positions = np.linspace(0, width - 1, len(gps_coords['x_ticks'])) - x_labels = [f"{lon:.5f}E" for lon, lat in gps_coords['x_ticks']] - ax.set_xticks(x_pixel_positions) - ax.set_xticklabels(x_labels, fontsize=7, rotation=30) - ax.set_xlabel('Longitude', fontsize=9, fontweight='bold') - - # Y-axis: latitude along left - y_pixel_positions = np.linspace(0, height - 1, len(gps_coords['y_ticks'])) - y_labels = [f"{lat:.5f}N" for lon, lat in gps_coords['y_ticks']] - ax.set_yticks(y_pixel_positions) - ax.set_yticklabels(y_labels, fontsize=7) - ax.set_ylabel('Latitude', fontsize=9, fontweight='bold') - else: - # Fallback to Lambert 93 coordinates - x_ticks_count = 5 - x_positions = np.linspace(0, width - 1, x_ticks_count) - x_labels = [f"{(min_x + xp * pixel_size_x)/1000:.1f}" for xp in x_positions] - ax.set_xticks(x_positions) - ax.set_xticklabels(x_labels, fontsize=8) - ax.set_xlabel('Est (km) - Lambert 93', fontsize=9, fontweight='bold') - - y_ticks_count = 5 - y_positions = np.linspace(0, height - 1, y_ticks_count) - y_labels = [f"{(max_y - yp * pixel_size_y)/1000:.1f}" for yp in y_positions] - ax.set_yticks(y_positions) - ax.set_yticklabels(y_labels, fontsize=8) - ax.set_ylabel('Nord (km) - Lambert 93', fontsize=9, fontweight='bold') - - # Style axes - ax.tick_params(axis='both', which='both', direction='out', length=3, - width=0.8, colors='black') - for spine in ax.spines.values(): - spine.set_visible(True) - spine.set_color('black') - spine.set_linewidth(0.8) - - # ---- North arrow ---- - from matplotlib.patches import Polygon as MplPolygon - from mpl_toolkits.axes_grid1.inset_locator import inset_axes - north_ax = inset_axes(ax, width="4%", height="7%", loc='upper right', - bbox_to_anchor=(-0.03, 0.12, 1, 1), bbox_transform=ax.transAxes) - north_ax.set_xlim(0, 1) - north_ax.set_ylim(0, 1) - north_ax.axis('off') - north_ax.plot([0.5, 0.5], [0.1, 0.65], color='black', linewidth=2.5, zorder=10) - north_ax.add_patch(MplPolygon([[0.5, 0.2], [0.35, 0.4], [0.5, 0.65], [0.65, 0.4]], - closed=True, facecolor='black', edgecolor='black', zorder=9)) - north_ax.text(0.5, 0.95, 'N', ha='center', va='top', - fontsize=13, fontweight='bold', color='black', zorder=11) - - # ---- Bottom cartouche: info + scale bar OUTSIDE the map ---- - info_ax = fig.add_subplot(gs[1]) - info_ax.axis('off') - - # Altitude range (skip for RGB images) - if is_rgb: - alt_min = 0 - alt_max = 0 - else: - alt_min = np.nanmin(valid_data) if len(valid_data) > 0 else 0 - alt_max = np.nanmax(valid_data) if len(valid_data) > 0 else 0 - extent_km_x = (max_x - min_x) / 1000 - extent_km_y = (max_y - min_y) / 1000 - - # Build info text (left side of bottom row) - if gps_coords: - nw_lat, nw_lon = gps_coords['NW'] - se_lat, se_lon = gps_coords['SE'] - info_text = ( - f"GPS: {nw_lat:.5f}N {nw_lon:.5f}E - {se_lat:.5f}N {se_lon:.5f}E | " - f"EPSG:2154 | Res: {self.resolution}m/px | " - f"Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km" - ) - else: - info_text = ( - f"EPSG:2154 | X: {min_x:.0f}-{max_x:.0f} Y: {min_y:.0f}-{max_y:.0f} | " - f"Res: {self.resolution}m/px | Emprise: {extent_km_x:.1f}x{extent_km_y:.1f}km" - ) - - info_ax.text(0.01, 0.5, info_text, - transform=info_ax.transAxes, fontsize=8.5, - verticalalignment='center', family='monospace', - bbox=dict(boxstyle='round,pad=0.3', facecolor='#f0f0f0', edgecolor='#aaaaaa', alpha=0.95)) - - # Scale bar: always 100m (fixed scale) - # 100m at 0.5m/px = 200 pixels - scale_m = 100 - scale_label = "100 m" - pixels_per_meter = 1.0 / pixel_size_x - scale_px = int(scale_m * pixels_per_meter) - # Convert pixel length to fraction of the axis width - # Scale bar as fraction of the info axis (which matches map width) - scale_bar_frac = scale_px / width - - # Draw scale bar as simple black bar with text - from matplotlib.patches import Rectangle as RectPatch - bar_left = 0.80 - bar_bottom = 0.15 - bar_width_frac = min(scale_bar_frac, 0.15) # Cap at 15% of width - bar_height = 0.35 - - info_ax.add_patch(RectPatch((bar_left, bar_bottom), bar_width_frac, bar_height, - facecolor='black', edgecolor='black', linewidth=0.5, - transform=info_ax.transAxes, clip_on=False)) - info_ax.text(bar_left + bar_width_frac / 2, bar_bottom + bar_height + 0.12, - scale_label, ha='center', va='bottom', fontsize=9, fontweight='bold', - transform=info_ax.transAxes) - - fig.patch.set_facecolor('white') - - # Save as PNG first (matplotlib doesn't support WebP well) - png_temp = self.vis_dir / f"{tif_file.stem}_temp.png" - plt.savefig(png_temp, dpi=150, bbox_inches='tight', pad_inches=0.15, - facecolor='white', format='png') - plt.close() - - # Convert PNG to lossless WebP using PIL - from PIL import Image as PILImage - img = PILImage.open(str(png_temp)) - img.save(str(jpg_file), format='WEBP', lossless=True) - png_temp.unlink() # Delete temp PNG - - # Delete the source TIFF file to save space - tif_file.unlink() - - return jpg_file - - except Exception as e: - print(f" Erreur conversion PNG: {e}") - import traceback - traceback.print_exc() - return None - - def generate_all_visualizations(self, dtm_file, basename): - """Generate all archaeological visualizations""" - print(f"\n Génération visualisations:") - - # Create per-file subdirectory - file_vis_dir = self.vis_dir / basename - file_vis_dir.mkdir(exist_ok=True) - # Temporarily override vis_dir for this file - original_vis_dir = self.vis_dir - self.vis_dir = file_vis_dir - - vis_results = {} - - # Core visualizations - vis_results['hillshade'] = self.generate_hillshade(dtm_file, basename) - vis_results['slope'] = self.generate_slope(dtm_file, basename) - vis_results['aspect'] = self.generate_aspect(dtm_file, basename) - vis_results['curvature'] = self.generate_curvature(dtm_file, basename) - vis_results['svf'] = self.generate_svf(dtm_file, basename) - vis_results['lrm'] = self.generate_lrm(dtm_file, basename) - vis_results['pos_open'] = self.generate_openness(dtm_file, basename, positive=True) - vis_results['neg_open'] = self.generate_openness(dtm_file, basename, positive=False) - - # Advanced visualizations - vis_results['mslrm'] = self.generate_mslrm(dtm_file, basename) - vis_results['tpi'] = self.generate_tpi(dtm_file, basename) - vis_results['depressions'] = self.generate_depressions(dtm_file, basename) - vis_results['sailore'] = self.generate_sailore(dtm_file, basename) - vis_results['roughness'] = self.generate_roughness(dtm_file, basename) - vis_results['anomalies'] = self.generate_anomalies(dtm_file, basename) - vis_results['wavelet'] = self.generate_wavelet(dtm_file, basename) - vis_results['texture'] = self.generate_texture(dtm_file, basename) - vis_results['flow'] = self.generate_flow(dtm_file, basename) - - # IGN overlays - vis_results['ortho'] = self.generate_ign_overlay( - dtm_file, basename, - layer='ORTHOIMAGERY.ORTHOPHOTOS', - title='Photographie Aerienne IGN', - legend_label='Orthophotographie\nImage aerienne', - description='Photographie aerienne IGN (Orthophoto)', - out_suffix='ortho' - ) - vis_results['topo'] = self.generate_ign_overlay( - dtm_file, basename, - layer='GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2', - title='Carte Topographique IGN', - legend_label='Carte IGN\nPlan topographique', - description='Carte topographique IGN (Plan IGN)', - out_suffix='topo' - ) - - # Convert to PNG - print(f"\n Conversion images PNG:") - for name, tif_file in vis_results.items(): - if tif_file: - jpg_file = self.tif_to_png(tif_file) - if jpg_file: - print(f" ✓ {jpg_file.name}") - - # Restore original vis_dir - self.vis_dir = original_vis_dir - - return vis_results - - # ============ STEP 5: PDF Report ============ - - def generate_pdf_report(self, basename): - """Generate A3 PDF report for a LiDAR file with all visualizations. - Page 1: Mise en situation (ortho + topo IGN side by side) - Pages 2+: Other visualizations (2 per page)""" - from matplotlib.backends.backend_pdf import PdfPages - - pdf_file = self.pdf_dir / f"{basename}_rapport.pdf" - print(f" → Génération rapport PDF A3: {pdf_file.name}") - - # Look for PNGs in per-file subdirectory first, then fallback to main dir - file_vis_dir = self.vis_dir / basename - if file_vis_dir.exists(): - png_files = sorted(file_vis_dir.glob("*.webp")) - else: - png_files = sorted(self.vis_dir.glob(f"{basename}_*.webp")) - if not png_files: - print(f" ✗ Aucune image PNG trouvée pour {basename}") - return None - - # Categorize files - situ_files = [] # Page 1: situation maps - analysis_files = [] # Pages 2+: analysis maps - - for f in png_files: - name = f.stem.lower() - if 'ortho' in name: - situ_files.insert(0, f) # Ortho first - elif 'topo' in name: - situ_files.append(f) # Topo second - else: - analysis_files.append(f) - - # Define display order for analysis maps - order = ['mslrm', 'svf', 'negative_openness', - 'positive_openness', 'sailore', 'depressions', 'hillshade_multi', - 'lrm', 'tpi', 'slope', 'curvature', 'aspect', - 'roughness', 'anomalies', 'wavelet', 'texture', 'flow'] - def sort_key(f): - name = f.stem.lower() - for i, key in enumerate(order): - if key in name: - return i - return len(order) - - analysis_files.sort(key=sort_key) - - # A3 landscape dimensions in inches (420mm x 297mm) - a3_w, a3_h = 16.54, 11.69 - - try: - with PdfPages(str(pdf_file)) as pdf: - # ============ PAGE 1: Mise en situation ============ - if situ_files: - fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white') - - n_situ = len(situ_files) - if n_situ == 2: - # Side by side - gs = fig.add_gridspec(1, 2, wspace=0.05, left=0.03, right=0.97, - top=0.92, bottom=0.06) - else: - # Single or more - gs = fig.add_gridspec(1, max(n_situ, 1), wspace=0.05, - left=0.03, right=0.97, top=0.92, bottom=0.06) - - fig.text(0.5, 0.97, f"Mise en situation - {basename}", - fontsize=20, fontweight='bold', ha='center', va='top') - - for i, f in enumerate(situ_files): - ax = fig.add_subplot(gs[0, i]) - img = plt.imread(str(f)) - ax.imshow(img) - ax.axis('off') - # Add title from filename - title = f.stem.replace(basename + '_', '').replace('_', ' ').title() - ax.set_title(title, fontsize=12, fontweight='bold', pad=5) - - pdf.savefig(fig, dpi=150) - plt.close(fig) - - # ============ PAGES 2+: Analysis maps (2 per page) ============ - # Group analysis files 2 per page - for page_start in range(0, len(analysis_files), 2): - page_files = analysis_files[page_start:page_start + 2] - - fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white') - - if len(page_files) == 2: - gs = fig.add_gridspec(1, 2, wspace=0.08, left=0.03, right=0.97, - top=0.93, bottom=0.05) - else: - gs = fig.add_gridspec(1, 1, left=0.05, right=0.95, - top=0.93, bottom=0.05) - - for i, f in enumerate(page_files): - ax = fig.add_subplot(gs[0, i]) - img = plt.imread(str(f)) - ax.imshow(img) - ax.axis('off') - # Title from filename - title = f.stem.replace(basename + '_', '').replace('_', ' ').title() - ax.set_title(title, fontsize=11, fontweight='bold', pad=3) - - # Page number - page_num = (page_start // 2) + 2 # Start from page 2 - fig.text(0.99, 0.01, f"Page {page_num}", fontsize=8, - ha='right', va='bottom', color='gray') - - pdf.savefig(fig, dpi=150) - plt.close(fig) - - print(f" ✓ Rapport PDF: {pdf_file.name}") - return pdf_file - - except Exception as e: - print(f" ✗ Erreur PDF: {e}") - import traceback - traceback.print_exc() - return None - - # ============ Complete Pipeline ============ - - def process_file(self, laz_file): - """Process a single LAZ file""" - basename = laz_file.stem - print(f"\n{'='*60}") - print(f" FICHIER : {basename}") - print(f"{'='*60}") - - # Step 1: Ground classification - print(f"\n[1/3] Classification du sol...") - las_file = self.classify_ground(laz_file) - if not las_file: - print(f" ✗ Échec classification") - return False - - # Step 2: Generate DTM - print(f"\n[2/3] Génération DTM...") - dtm_file = self.create_dtm_fast(las_file, basename) - if not dtm_file: - print(f" ✗ Échec DTM") - return False - - # Step 3: Visualizations - print(f"\n[3/4] Visualisations archéologiques...") - vis = self.generate_all_visualizations(dtm_file, basename) - - # Step 4: PDF report - print(f"\n[4/4] Rapport PDF A3...") - self.generate_pdf_report(basename) - - print(f"\n✓ {basename} traité avec succès !") - return True - - def process_all(self): - """Process all LAZ files in input directory""" - files = self.find_laz_files() - - if not files: - print("Aucun fichier LAZ/LAS trouvé !") - return - - print(f"\n{'='*60}") - print(f" PIPELINE ARCHÉOLOGIQUE LiDAR") - print(f"{'='*60}") - - # Check tools - print(f"\nVérification des outils...") - if not self.check_tools(): - print("Erreur: Certains outils ne sont pas disponibles") - return - - results = {} - - if self.workers > 1 and len(files) > 1: - # Process multiple files in parallel - print(f"\n Traitement parallèle avec {self.workers} workers...") - with ProcessPoolExecutor(max_workers=self.workers) as executor: - future_to_file = { - executor.submit(self._process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), self.resolution): laz_file - for laz_file in files - } - for future in as_completed(future_to_file): - laz_file = future_to_file[future] - try: - success = future.result() - results[laz_file.name] = success - status = "✓" if success else "✗" - print(f" {status} {laz_file.name}") - except Exception as e: - print(f"\n✗ Erreur traitement {laz_file.name}: {e}") - results[laz_file.name] = False - else: - # Sequential processing - for laz_file in files: - try: - results[laz_file.name] = self.process_file(laz_file) - except Exception as e: - print(f"\n✗ Erreur traitement: {e}") - import traceback - traceback.print_exc() - results[laz_file.name] = False - - # Summary - print(f"\n{'='*60}") - print(f" RÉSUMÉ") - print(f"{'='*60}") - success_count = sum(results.values()) - print(f"✓ {success_count}/{len(results)} fichiers traités avec succès") - - print(f"\nRésultats dans: {self.output_dir}") - print(f" • DTM: {self.dtm_dir}") - print(f" • Visualisations PNG: {self.vis_dir}") - print(f" • Rapports PDF: {self.pdf_dir}") - - # Clean up temporary files to save space - print(f"\nNettoyage des fichiers temporaires...") - try: - import shutil - if self.temp_dir.exists(): - shutil.rmtree(self.temp_dir) - print(f" ✓ Fichiers temporaires supprimés ({self.temp_dir})") - except Exception as e: - print(f" Note: Impossible de supprimer les fichiers temporaires") - - @staticmethod - def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution): - """Standalone function for multiprocessing - creates its own pipeline instance.""" - pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1) - laz_file = Path(laz_file_str) - return pipeline.process_file(laz_file) - - -def main(): - import argparse - - parser = argparse.ArgumentParser( - description="Pipeline LiDAR pour détection archéologique (Python pur)" - ) - parser.add_argument( - "input", - help="Dossier contenant les fichiers LAZ/LAS" - ) - parser.add_argument( - "-o", "--output", - default="/data/output", - help="Dossier de sortie" - ) - parser.add_argument( - "-r", "--resolution", - type=float, - default=0.5, - help="Résolution en mètres" - ) - parser.add_argument( - "-w", "--workers", - type=int, - default=1, - help="Nombre de workers pour traitement parallèle (défaut: 1)" - ) - - args = parser.parse_args() - - try: - pipeline = LidarArchaeoPipeline( - input_dir=args.input, - output_dir=args.output, - resolution=args.resolution, - workers=args.workers - ) - pipeline.process_all() - except Exception as e: - print(f"Erreur: {e}") - sys.exit(1) - +from lidar_pipeline.cli import main if __name__ == "__main__": - main() + main() \ No newline at end of file diff --git a/run.sh b/run.sh index aceec6e..018a9cf 100755 --- a/run.sh +++ b/run.sh @@ -1,11 +1,16 @@ #!/bin/bash # Pipeline LiDAR Archéologique # Utilisation: ./run.sh [options] +# # Options: # -r RESOLUTION Résolution en m/px (défaut: 0.5) # -w WORKERS Nombre de workers parallèles (défaut: 1) # -g Activer l'accélération GPU -# -h Afficher l'aide +# -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 +# --file NOM Traiter un seul fichier LAZ (pour tests) +# -h Afficher l'aide complète set -e @@ -16,33 +21,59 @@ IMAGE_NAME="lidar-lidar" RESOLUTION=0.5 WORKERS=1 GPU_FLAG="" +VERBOSE_FLAG="" +FORCE_FLAG="" +FILE_FILTER="" -while getopts "r:w:gh" opt; do +while getopts "r:w:gvf-:" opt; do case $opt in r) RESOLUTION="$OPTARG" ;; w) WORKERS="$OPTARG" ;; g) GPU_FLAG="--gpus all" ;; + v) VERBOSE_FLAG="-v" ;; + f) FORCE_FLAG="--force" ;; + -) + case "${OPTARG}" in + debug) VERBOSE_FLAG="--debug" ;; + force) FORCE_FLAG="--force" ;; + file) ;; + *) echo "Option invalide: --${OPTARG}" >&2; exit 1 ;; + esac + ;; h) echo "Pipeline LiDAR Archéologique" echo "" - echo "Usage: $0 [-r RESOLUTION] [-w WORKERS] [-g]" + echo "Usage: $0 [options]" echo "" echo " -r RESOLUTION Résolution en m/px (défaut: 0.5)" 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)" + echo " --debug Mode debug (détails internes fichier:ligne)" + echo " -f / --force Régénérer tous les fichiers même si les WebP existent" + echo " --file NOM Traiter un seul fichier LAZ (pour tests)" echo " -h Afficher cette aide" echo "" echo "Exemples:" echo " $0 # Traitement standard" echo " $0 -g # Avec accélération GPU" echo " $0 -g -w 4 # GPU + 4 workers parallèles" - echo " $0 -r 0.2 -g # Haute résolution + GPU" + echo " $0 -g -v # GPU + mode verbeux" + echo " $0 -r 0.2 -g --debug # Haute résolution + GPU + debug" + echo " $0 -f # Forcer la régénération de tous les fichiers" + echo " $0 --file 6881 # Traiter un seul fichier (pour tests)" exit 0 ;; *) echo "Option invalide. Utilisez -h pour l'aide." >&2; exit 1 ;; esac done +# Handle --file option separately (not easily done in getopts) +if [[ " $@ " == *" --file "* ]]; then + # Extract the argument after --file + FILE_FILTER=$(echo "$@" | sed -n 's/.*--file *\([^ ]*\).*/\1/p') +fi + # Build l'image si elle n'existe pas if ! docker image inspect "$IMAGE_NAME" >/dev/null 2>&1; then echo "Build de l'image Docker..." @@ -59,14 +90,22 @@ echo "============================================" echo " Résolution : ${RESOLUTION}m/px" echo " Workers : ${WORKERS}" echo " GPU : $([ -n "$GPU_FLAG" ] && echo 'OUI' || echo 'non')" +echo " Verbeux : $([ -n "$VERBOSE_FLAG" ] && echo 'OUI' || echo 'non')" +echo " Force : $([ -n "$FORCE_FLAG" ] && echo 'OUI' || echo 'non')" +if [ -n "$FILE_FILTER" ]; then + echo " Fichier : ${FILE_FILTER}" +fi echo "============================================" +CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG" +if [ -n "$FILE_FILTER" ]; then + CMD_ARGS="$CMD_ARGS --file $FILE_FILTER" +fi + docker run --rm $GPU_FLAG \ --user 1000:1000 \ -v "${INPUT_DIR}:/data/input:ro" \ -v "${OUTPUT_DIR}:/data/output" \ "$IMAGE_NAME" \ - python3 /usr/local/bin/process_lidar.py /data/input \ - -o /data/output \ - -r "$RESOLUTION" \ - -w "$WORKERS" \ No newline at end of file + python3 -m lidar_pipeline /data/input \ + $CMD_ARGS \ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..7e8ff94 --- /dev/null +++ b/setup.py @@ -0,0 +1,15 @@ +"""Setup for lidar_pipeline package.""" +from setuptools import setup, find_packages + +setup( + name='lidar_pipeline', + version='2.0.0', + description='Pipeline LiDAR pour détection archéologique', + packages=find_packages(), + python_requires='>=3.8', + entry_points={ + 'console_scripts': [ + 'lidar-pipeline=lidar_pipeline.cli:main', + ], + }, +) \ No newline at end of file