diff --git a/Dockerfile b/Dockerfile index 024b195..2406fbc 100644 --- a/Dockerfile +++ b/Dockerfile @@ -26,15 +26,20 @@ RUN pip3 install --no-cache-dir \ 'laspy[laspy]' \ scikit-image \ scikit-learn \ + scipy \ tqdm # Copy scripts COPY process_lidar.py /usr/local/bin/ RUN chmod +x /usr/local/bin/process_lidar.py -# Create directories with correct permissions -RUN mkdir -p /data/output /data/input && \ - chmod 777 /data /data/output /data/input +# Create user with uid/gid 1000:1000 and run as that user +RUN groupadd -g 1000 lidar && \ + useradd -u 1000 -g lidar -m lidar && \ + mkdir -p /data/output /data/input && \ + chown -R lidar:lidar /data /data/output /data/input + +USER lidar VOLUME ["/data"] diff --git a/README.md b/README.md index fe10479..27622b4 100644 --- a/README.md +++ b/README.md @@ -2,180 +2,185 @@ Workflow automatisé pour générer des visualisations exploitables à partir de données LiDAR pour la détection de structures archéologiques. -## 🎯 Ce que détecte ce pipeline +## Visualisations (21 par fichier) -| Visualisation | Utilité archéologique | -|--------------|----------------------| -| **Hillshade multidirectionnel** | Murs, terrasses, structures linéaires, routes | -| **Slope (Pente)** | Murs de soutènement, talus, changements brusques | -| **Sky-View Factor** | ⭐ **Top** - structures, tumulus, fondations | -| **Local Relief Model** | Micro-reliefs, fossés, levées de terrain | -| **Positive Openness** | Élévations, tumulus, bâtiments | -| **Negative Openness** | ⭐ **Cavités, fossés, souterrains** | +### Visualisations principales +| # | Visualisation | Utilité archéologique | +|---|--------------|----------------------| +| 1 | **Hillshade multidirectionnel** | Murs, terrasses, structures linéaires, routes | +| 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) | -## 📦 Installation Docker +### 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 | +| 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 | + +### Cartes de référence IGN +| # | Visualisation | Source | +|---|--------------|--------| +| 20 | **Photographie aérienne IGN** | Orthophotographie WMTS | +| 21 | **Carte topographique IGN** | Plan IGN V2 WMTS | + +## Installation Docker ```bash -# Clone ou copiez les fichiers dans un dossier cd /votre/dossier/lidar - -# Créez le dossier pour vos fichiers LAZ mkdir -p input # Copiez vos fichiers .laz dans input/ cp /chemin/vos/fichiers/*.laz input/ # Build l'image Docker -docker-compose build - -# Ou avec docker build docker build -t lidar-archeo . ``` -## 🚀 Utilisation - -### Méthode 1: docker-compose (recommandé) +## Utilisation +### Traitement standard (1 CPU) ```bash -# Traitement avec résolution 0.5m -docker-compose up - -# Résolution plus fine (plus lent) -docker-compose run -e RESOLUTION=0.2 lidar process_lidar.py /data/input -o /data/output -r 0.2 -``` - -### Méthode 2: docker run - -```bash -# Basic docker run --rm \ -v $(pwd)/input:/data/input:ro \ -v $(pwd)/output:/data/output \ lidar-archeo +``` -# Avec résolution personnalisée +### 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 -r 0.3 + process_lidar.py /data/input -o /data/output -w 4 +``` -# Plus de mémoire pour fichiers volumineux +### Résolution personnalisée +```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 + +# Résolution standard (0.5m) - recommandée archéologie docker run --rm \ - --memory=16g \ - --memory-swap=16g \ -v $(pwd)/input:/data/input:ro \ -v $(pwd)/output:/data/output \ lidar-archeo ``` -## 📁 Structure des dossiers +### Combinaison résolution + multi-CPU +```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 +``` + +## Structure des dossiers ``` . -├── input/ # Vos fichiers .laz (monté en read-only) -├── output/ # Résultats générés -│ ├── DTM/ # Modèles numériques d'élévation -│ ├── visualisations/# Images PNG prêtes à l'emploi -│ └── rapports/ # Vues synthétiques +├── input/ # Vos 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 +│ │ ├── fichier_6881/ # Un sous-dossier par fichier LAZ +│ │ │ ├── ..._hillshade_multi.png +│ │ │ ├── ..._svf.png +│ │ │ ├── ..._mslrm.png +│ │ │ └── ... (21 visualisations) +│ │ └── fichier_6882/ +│ │ └── ... +│ └── rapports/ # Rapports PDF A3 par fichier +│ ├── fichier_6881_rapport.pdf +│ └── fichier_6882_rapport.pdf ├── Dockerfile ├── docker-compose.yml └── process_lidar.py ``` -## 📊 Sortie générée +## Paramètres -Après traitement, dans `output/` : +| Paramètre | Option | Défaut | Description | +|-----------|--------|--------|-------------| +| 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 | -``` -output/ -├── DTM/ -│ └── fichier_dtm.tif -├── visualisations/ -│ ├── fichier_hillshade_multi.png -│ ├── fichier_svf.png # ⭐ Sky-View Factor -│ ├── fichier_lrm.png # Local Relief Model -│ ├── file_positive_openness.png # Élévations -│ ├── file_negative_openness.png # ⭐ Cavités -│ └── fichier_slope.png -└── rapports/ - └── fichier_overview.png # Vue 2x3 synthétique -``` +### 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 -## 👁️ Interprétation +## Interprétation archéologique ### Pour détecter les cavités et souterrains -Regardez dans cet ordre : -1. **Negative Openness** - Zones bleues foncées = creux -2. **Local Relief Model** - Zones bleues = dépressions -3. **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. **Sky-View Factor** - Structures géométriques claires -2. **Positive Openness** - Zones orangées/rouges = élévations -3. **Hillshade** - Lignes droites, rectangles, angles +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 élévations de terrain -1. **Slope** - Changements de pente brusques -2. **Positive Openness** - Zones en hauteur -3. **LRM** - Zones rouges = relief local +### 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 -## ⚙️ Paramètres +### 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 -### Résolution (-r) -- `0.2` - Très fine, bâtiments individuels (lent) -- `0.5` - Recommandé archéologie (équilibre) -- `1.0` - Rapide, grandes structures - -### Ressources Docker -Modifiez dans `docker-compose.yml` : -```yaml -deploy: - resources: - limits: - cpus: '4' # CPU cores - memory: 8G # RAM -``` - -## 🔧 Dépannage +## Dépannage ```bash -# Vérifier que Docker fonctionne +# Vérifier Docker docker --version -docker-compose --version - -# Voir les logs en temps réel -docker-compose up -f # Shell dans le conteneur -docker-compose run lidar bash +docker run --rm -it \ + -v $(pwd)/input:/data/input \ + -v $(pwd)/output:/data/output \ + --entrypoint bash lidar-archeo -# Nettoyer tout -docker-compose down -v +# Reconstruire l'image +docker build --no-cache -t lidar-archeo . + +# Nettoyer docker system prune -a ``` ### Erreur mémoire -```bash -# Augmenter la limite mémoire dans docker-compose.yml -# Ou utiliser --memory=16g avec docker run -``` +Augmenter la mémoire Docker à 16Go+ pour les gros fichiers LiDAR HD. -### Erreur PDAL/Whitebox -```bash -# Recréer l'image -docker-compose build --no-cache -``` - -## 📝 Prochaines améliorations - -- [ ] Classification automatique des structures détectées -- [ ] Export GeoTIFF géoréférencés pour QGIS -- [ ] Interface web pour exploration interactive -- [ ] Support multi-fichiers avec mosaïque - -## 📚 Ressources +## Ressources - [PDAL Documentation](https://pdal.io/) -- [WhiteboxTools](https://www.whiteboxgeo.com/manual/wbt_book/) -- [Archéologie LiDAR](https://archaeologydataservice.ac.uk/) +- [LiDAR Archéologie - Méthodes avancées](https://archaeologydataservice.ac.uk/) +- [Relief Visualization Toolbox](https://rvtpy.readthedocs.io/) \ No newline at end of file diff --git a/process_lidar.py b/process_lidar.py index 6fd27a7..cbb59a6 100755 --- a/process_lidar.py +++ b/process_lidar.py @@ -8,6 +8,7 @@ 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 @@ -26,6 +27,12 @@ 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 + rcParams['figure.dpi'] = 150 rcParams['savefig.dpi'] = 300 rcParams['font.size'] = 10 @@ -34,10 +41,11 @@ rcParams['font.size'] = 10 class LidarArchaeoPipeline: """Pipeline Python pur pour analyse archéologique LiDAR""" - def __init__(self, input_dir, output_dir, resolution=0.5): + 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(): @@ -48,14 +56,16 @@ class LidarArchaeoPipeline: self.dtm_dir = self.output_dir / "DTM" self.vis_dir = self.output_dir / "visualisations" + self.pdf_dir = self.output_dir / "rapports" - for d in [self.dtm_dir, self.vis_dir]: + 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""" @@ -233,38 +243,47 @@ class LidarArchaeoPipeline: # Fill gaps using interpolation print(f" Remplissage des zones vides...") - from scipy.ndimage import distance_transform_edt # Create mask of empty cells mask = np.isnan(dtm) if np.any(mask): - # Use inpainting to fill gaps + 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 - # Get coordinates of valid points y_coords, x_coords = np.where(~mask) z_values = dtm[~mask] - # Create interpolation function - if len(y_coords) > 100: # Only interpolate if enough points - # Create grid for interpolation - y_all, x_all = np.mgrid[0:dtm.shape[0], 0:dtm.shape[1]] - - # Use nearest neighbor interpolation for speed + if len(y_coords) > 100: + # Nearest neighbor interpolation (fast, memory-efficient) interp = interpolate.NearestNDInterpolator( np.column_stack((y_coords, x_coords)), z_values ) - # Fill missing values - dtm_filled = interp(y_all, x_all) + # Only interpolate where needed (memory-efficient) + y_missing, x_missing = np.where(mask) + dtm[y_missing, x_missing] = interp(y_missing, x_missing) - # Apply slight smoothing to blend filled areas - dtm_filled = gaussian_filter(dtm_filled, sigma=0.5) + # 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) - # Replace NaN values with filled values - dtm[mask] = dtm_filled[mask] + # 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" @@ -321,7 +340,7 @@ class LidarArchaeoPipeline: # Calculate slope and aspect slope = np.arctan(np.sqrt(dx**2 + dy**2)) - aspect = np.arctan2(-dy, dx) + aspect = np.arctan2(dy, dx) # Hillshade formula hs = np.sin(alt_rad) * np.sin(slope) + \ @@ -351,6 +370,56 @@ class LidarArchaeoPipeline: 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)...") @@ -401,7 +470,7 @@ class LidarArchaeoPipeline: # Calculate aspect (direction of slope) dy, dx = np.gradient(dem) - aspect = np.arctan2(-dy, dx) * 180 / np.pi + aspect = np.arctan2(dy, dx) * 180 / np.pi aspect = np.mod(aspect, 360) # Convert to 0-360 # Save @@ -482,7 +551,7 @@ class LidarArchaeoPipeline: # Calculate slope and aspect slope = np.arctan(np.sqrt(dx**2 + dy**2)) - aspect = np.arctan2(-dy, dx) + aspect = np.arctan2(dy, dx) # Solar irradiance (morning sun - azimuth 90, altitude 30) az_rad = np.radians(90) @@ -555,8 +624,11 @@ class LidarArchaeoPipeline: return None def generate_svf(self, dem_file, basename): - """Sky-View Factor approximation""" - print(f" → Sky-View Factor...") + """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.""" + print(f" → Sky-View Factor (ray-tracing)...") output = self.vis_dir / f"{basename}_svf.tif" @@ -566,12 +638,59 @@ class LidarArchaeoPipeline: transform = src.transform crs = src.crs - # Simplified SVF using relief inverted - # Normalize DEM - dem_norm = (dem - np.nanmin(dem)) / (np.nanmax(dem) - np.nanmin(dem)) + rows, cols = dem.shape + res = self.resolution # meters per pixel - # Apply inverted relief (approximation of SVF) - svf = 1 - dem_norm + # 16 azimuth directions (0°, 22.5°, 45°, ... 337.5°) + n_dirs = 16 + angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) + dx = np.cos(angles) + dy = np.sin(angles) + + # Maximum search distance (in pixels) + max_dist = int(50 / res) # 50m search radius + + # Pad DEM with NaN to handle boundaries + padded = np.pad(dem.astype(np.float64), max_dist, mode='constant', + constant_values=np.nan) + + svf = np.zeros_like(dem, dtype=np.float64) + + for d_idx in range(n_dirs): + # Direction vector + ddx, ddy = dx[d_idx], dy[d_idx] + + # For each step along the ray, compute the horizon angle + horizon = np.zeros_like(dem, dtype=np.float64) + + for step in range(1, max_dist + 1): + # Offset in pixels (rounded to nearest integer) + px = int(round(ddx * step)) + py = int(round(ddy * step)) + + # Distance in meters + dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) + if dist_m < res * 0.5: + continue + + # Elevation difference relative to center pixel + # Source is at (max_dist + row, max_dist + col) in padded array + # Target is at (max_dist + row + py, max_dist + col + px) + elev_diff = padded[max_dist + py:max_dist + py + rows, + max_dist + px:max_dist + px + cols] - dem + + # Horizon angle = arctan(elev_diff / dist) + angle = np.arctan2(elev_diff, dist_m) + + # Keep maximum horizon angle + horizon = np.where(np.isnan(angle), horizon, + np.maximum(horizon, np.nan_to_num(angle, nan=0))) + + # SVF contribution: cos²(zenith) where zenith = pi/2 - horizon + # Higher horizon = less sky visible + svf += np.cos(np.pi / 2 - horizon) ** 2 + + svf /= n_dirs # Average over all directions # Save with rasterio.open( @@ -594,9 +713,13 @@ class LidarArchaeoPipeline: return None def generate_openness(self, dem_file, basename, positive=True): - """Positive/Negative Openness approximation""" + """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.""" name = "positive_openness" if positive else "negative_openness" - print(f" → {name.replace('_', ' ').title()}...") + print(f" → {name.replace('_', ' ').title()} (ray-tracing)...") output = self.vis_dir / f"{basename}_{name}.tif" @@ -606,17 +729,1002 @@ class LidarArchaeoPipeline: transform = src.transform crs = src.crs - if positive: - # Positive openness: high points - from scipy.ndimage import maximum_filter - openness = maximum_filter(dem, size=int(10/self.resolution)) - else: - # Negative openness: low points (cavities!) - from scipy.ndimage import minimum_filter - openness = minimum_filter(dem, size=int(10/self.resolution)) + rows, cols = dem.shape + res = self.resolution - # Calculate difference - result = openness - dem if positive else dem - openness + # 8 cardinal directions + 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) # 50m search radius + + # Pad DEM with NaN to handle boundaries + padded = np.pad(dem.astype(np.float64), max_dist, mode='constant', + constant_values=np.nan) + + openness_sum = np.zeros_like(dem, dtype=np.float64) + + for d_idx in range(n_dirs): + ddx, ddy = dx[d_idx], dy[d_idx] + + # For positive openness: find max upward angle (zenith) + # For negative openness: find max downward angle (nadir) + max_angle = np.zeros_like(dem, dtype=np.float64) + + 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: + # Positive openness: angle from vertical to terrain above + # Only consider points higher than center (elev_diff > 0) + angle = np.arctan2(np.maximum(elev_diff, 0), dist_m) + else: + # Negative openness: angle from vertical to terrain below + # Only consider points lower than center (elev_diff < 0) + angle = np.arctan2(np.maximum(-elev_diff, 0), dist_m) + + max_angle = np.where(np.isnan(angle), max_angle, + np.maximum(max_angle, np.nan_to_num(angle, nan=0))) + + openness_sum += max_angle + + # Average across all directions (convert to degrees) + openness_result = np.degrees(openness_sum / n_dirs) + + # Save + 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).""" + print(f" → Multi-Scale Relief Model (MSRM)...") + + output = self.vis_dir / f"{basename}_mslrm.tif" + + try: + with rasterio.open(dem_file) as src: + dem = src.read(1) + transform = src.transform + crs = src.crs + + from scipy.ndimage import gaussian_filter + + # Compute LRM at multiple scales + sigmas = [5, 10, 25, 50, 100] + lrm_stack = [] + + for sigma in sigmas: + sigma_px = sigma / self.resolution + local_mean = gaussian_filter(dem, sigma=sigma_px) + lrm = dem - local_mean + # Normalize each scale + lrm_norm = lrm / max(np.nanstd(lrm), 0.01) + lrm_stack.append(lrm_norm) + + # Combine: RMS of normalized LRM at all scales + mslrm = np.sqrt(np.mean(np.array(lrm_stack) ** 2, axis=0)) + + with rasterio.open( + output, + 'w', + driver='GTiff', + height=mslrm.shape[0], + width=mslrm.shape[1], + count=1, + dtype='float32', + crs=crs, + transform=transform, + compress='lzw' + ) as dst: + dst.write(mslrm.astype('float32'), 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.""" + print(f" → TPI multi-echelle...") + + output = self.vis_dir / f"{basename}_tpi.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 + + # Fine-scale TPI (5m radius) + fine_size = int(5 / self.resolution) + if fine_size % 2 == 0: + fine_size += 1 + tpi_fine = dem - uniform_filter(dem, size=fine_size) + + # Broad-scale TPI (100m radius) + broad_size = int(100 / self.resolution) + if broad_size % 2 == 0: + broad_size += 1 + tpi_broad = dem - 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 + fine_std = max(np.nanstd(tpi_fine), 0.01) + broad_std = max(np.nanstd(tpi_broad), 0.01) + tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std) + + with rasterio.open( + output, + 'w', + driver='GTiff', + height=tpi_combined.shape[0], + width=tpi_combined.shape[1], + count=1, + dtype='float32', + crs=crs, + transform=transform, + compress='lzw' + ) as dst: + dst.write(tpi_combined.astype('float32'), 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.""" + print(f" → SAILORE (LRM adaptatif)...") + + output = self.vis_dir / f"{basename}_sailore.tif" + + try: + with rasterio.open(dem_file) as src: + dem = src.read(1) + transform = src.transform + crs = src.crs + + from scipy.ndimage import gaussian_filter, uniform_filter + + # Compute slope for adaptive kernel sizing + gy, gx = np.gradient(dem, self.resolution) + slope = np.arctan(np.sqrt(gx**2 + gy**2)) + slope_deg = np.degrees(slope) + + # Adaptive sigma: flat terrain (low slope) = large kernel, steep = small + # slope 0° -> sigma_max (25m), slope 30° -> sigma_min (2m) + sigma_min = 2.0 / self.resolution + sigma_max = 25.0 / self.resolution + # Normalize slope to 0-1 range + slope_norm = np.clip(slope_deg / 30.0, 0, 1) + # Invert: flat areas get high sigma, steep areas get low sigma + adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min) + + # Compute local mean with adaptive Gaussian (approximated by blending) + # For efficiency, compute at 3 fixed scales and blend based on slope + lrm_fine = dem - gaussian_filter(dem, sigma=sigma_min) + lrm_medium = dem - gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2) + lrm_coarse = dem - gaussian_filter(dem, sigma=sigma_max) + + # Blend based on slope: steep -> fine, flat -> coarse + w_fine = slope_norm + w_medium = 1 - 2 * np.abs(slope_norm - 0.5) + w_coarse = 1 - slope_norm + # Normalize weights + 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 + + with rasterio.open( + output, + 'w', + driver='GTiff', + height=sailore.shape[0], + width=sailore.shape[1], + count=1, + dtype='float32', + crs=crs, + transform=transform, + compress='lzw' + ) as dst: + dst.write(sailore.astype('float32'), 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.""" + print(f" → Ondelette Mexican Hat multi-echelle...") + + output = self.vis_dir / f"{basename}_wavelet.tif" + + try: + with rasterio.open(dem_file) as src: + dem = src.read(1) + transform = src.transform + crs = src.crs + + # Mexican Hat (Ricker) wavelet at multiple scales + # Uses scipy.ndimage.gaussian_laplace as 2D Mexican Hat approximation + scales = [2, 5, 10, 20, 50] # meters + wavelet_stack = [] + + for scale_m in scales: + # Create 1D Ricker wavelet and apply as 2D separable filter + sigma_px = scale_m / self.resolution + # 2D Mexican Hat = Laplacian of Gaussian + from scipy.ndimage import gaussian_laplace + response = -gaussian_laplace(dem.astype(np.float64), sigma=sigma_px) + # Normalize by scale + response /= max(np.nanstd(response), 0.01) + wavelet_stack.append(response) + + # Combine: RMS of all scales + combined = np.sqrt(np.mean(np.array(wavelet_stack) ** 2, axis=0)) + + 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 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( @@ -635,165 +1743,596 @@ class LidarArchaeoPipeline: return output except Exception as e: - print(f" ✗ Erreur openness: {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 JPEG with explicit legend""" + """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}.jpg" + jpg_file = self.vis_dir / f"{tif_file.stem}.png" try: with rasterio.open(tif_file) as src: - data = src.read(1) - nodata = src.nodata + # 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 nodata is not None: + 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) - # Get valid data for normalization - valid_data = data.compressed() if hasattr(data, 'compressed') else data.flatten() - valid_data = valid_data[~np.isnan(valid_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 'hillshade' in str(tif_file): + 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 (0-1)" - description = "Ombres → Structures, murs, terrasses visibles" + 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 (Slope)" - legend_label = f"Pente (°)\nMin: {vmin:.1f}° | Max: {vmax:.1f}°" - description = "Orange/Clair = Forte pente (murs, talus) | Foncé = Faible pente" + 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' # Cyclic colormap for directions - # Aspect is 0-360, normalize to 0-1 + cmap = 'hsv' vmin, vmax = 0, 360 data = np.clip((data - vmin) / (vmax - vmin), 0, 1) - title = "Aspect (Orientation)" - legend_label = "Orientation (° du Nord)\nN=0°, E=90°, S=180°, O=270°" - description = "Couleur = Direction de la pente (utile pour orientation bâtiments)" + 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' # Diverging for positive/negative curvature + 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 (Curvature)" - legend_label = f"Courbure\nRouge = Convexe (bosse)\nBleu = Concave (creux)" - description = "⭐ Excellent pour fossés, levées, terrasses, talus" - elif 'solar' in str(tif_file): - cmap = 'YlOrBr' # Sun-like colormap - vmin, vmax = 0, 1 - data = np.clip(data, vmin, vmax) - title = "Éclairage Solaire (Matin)" - legend_label = "Irradiance Solaire\nFoncé = Ombre | Clair = Éclairé" - description = "Simulation soleil matin - révèle structures orientées Est" + 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" - legend_label = "Ouverture vers le ciel\nFoncé = Cuvette | Clair = Sommet" - description = "Excellent pour structures, tumulus, fondations" + 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 = "Local Relief Model (LRM)" - legend_label = f"Relief local (m)\nRouge = +{vmax:.1f}m (élévation)\nBleu = {vmin:.1f}m (dépression)" - description = "Rouge = Surélévation | Bleu = Creux, fossé" + 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 = "Positive Openness" - legend_label = f"Ouverture positive (m)\nMin: {vmin:.1f}m | Max: {vmax:.1f}m" - description = "Orange = Zones élevées (tumulus, collines)" + 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 = "Negative Openness" - legend_label = f"Ouverture négative (m)\nMin: {vmin:.1f}m | Max: {vmax:.1f}m" - description = "Vert/Bleu = Zones basses (cavités, fossés, souterrains)" + 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 normalisée" + legend_label = "Altitude normalisee" description = "" - # Create figure with white background for JPEG - fig, ax = plt.subplots(figsize=(18, 12), facecolor='white') + # ---- 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) - # Display data with north at the top - im = ax.imshow(data, cmap=cmap, aspect='equal', origin='upper') + # 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') - # Enhanced colorbar with explicit values - cbar = plt.colorbar(im, ax=ax, pad=0.03, shrink=0.75, aspect=30) - cbar.ax.tick_params(labelsize=10, width=1.5) - cbar.outline.set_linewidth(1.5) + # Title above the map + ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10) - # Add colorbar label - cbar.set_label(legend_label, fontsize=11, fontweight='bold') + # 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') - # Enhanced title with description - full_title = f"{title}\n{description}" - ax.set_title(full_title, fontsize=16, fontweight='bold', pad=15) + # ---- 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') - # Add coordinates info - ax.text(0.02, 0.98, f"Résolution: {self.resolution}m/px", - transform=ax.transAxes, fontsize=9, - verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.8)) + # 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') - ax.axis('off') + 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') - # Add vector north arrow above the colorbar on the right side (black lines) - from matplotlib.patches import FancyArrow, Polygon + # 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 - - # Create inset axes positioned above the colorbar on the right - # The colorbar is typically at the right, so we place north arrow above it - north_ax = inset_axes(ax, width="5%", height="8%", loc='upper right', - bbox_to_anchor=(-0.02, 0.15, 1, 1), bbox_transform=ax.transAxes) - + 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') - - # Draw vector arrow pointing north (black lines) - # Main arrow shaft - north_ax.plot([0.5, 0.5], [0.15, 0.65], color='black', linewidth=2.5, zorder=10) - # Arrow head (triangle outline) - arrow_head_outline = [[0.5, 0.25], [0.35, 0.45], [0.5, 0.65], [0.65, 0.45]] - for i in range(len(arrow_head_outline) - 1): - north_ax.plot([arrow_head_outline[i][0], arrow_head_outline[i+1][0]], - [arrow_head_outline[i][1], arrow_head_outline[i+1][1]], - color='black', linewidth=2.5, zorder=10) - # Fill the arrow head - north_ax.add_patch(Polygon([[0.5, 0.25], [0.35, 0.45], [0.5, 0.65], [0.65, 0.45]], + 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) - # Add "N" label above the arrow - north_ax.text(0.5, 0.92, 'N', ha='center', va='top', - fontsize=14, 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') - plt.tight_layout() - # Save as JPEG - plt.savefig(jpg_file, dpi=150, bbox_inches='tight', pad_inches=0.1, facecolor='white', format='jpg') + plt.savefig(jpg_file, dpi=150, bbox_inches='tight', pad_inches=0.15, + facecolor='white', format='png') plt.close() # Delete the source TIFF file to save space @@ -802,7 +2341,7 @@ class LidarArchaeoPipeline: return jpg_file except Exception as e: - print(f" ✗ Erreur conversion JPEG: {e}") + print(f" Erreur conversion PNG: {e}") import traceback traceback.print_exc() return None @@ -811,28 +2350,189 @@ class LidarArchaeoPipeline: """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 = {} - # Generate rasters (existing + new) + # 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['solar'] = self.generate_solar(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) - # Convert to JPEG - print(f"\n Conversion images JPEG:") + # 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(): - jpg_file = self.tif_to_png(tif_file) - if jpg_file: - print(f" ✓ {jpg_file.name}") + 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("*.png")) + else: + png_files = sorted(self.vis_dir.glob(f"{basename}_*.png")) + 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): @@ -857,9 +2557,13 @@ class LidarArchaeoPipeline: return False # Step 3: Visualizations - print(f"\n[3/3] Visualisations archéologiques...") + 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 @@ -882,14 +2586,35 @@ class LidarArchaeoPipeline: return results = {} - 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 + + 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}") @@ -900,7 +2625,8 @@ class LidarArchaeoPipeline: print(f"\nRésultats dans: {self.output_dir}") print(f" • DTM: {self.dtm_dir}") - print(f" • Visualisations JPEG: {self.vis_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...") @@ -912,6 +2638,13 @@ class LidarArchaeoPipeline: 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 @@ -934,6 +2667,12 @@ def main(): 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() @@ -941,7 +2680,8 @@ def main(): pipeline = LidarArchaeoPipeline( input_dir=args.input, output_dir=args.output, - resolution=args.resolution + resolution=args.resolution, + workers=args.workers ) pipeline.process_all() except Exception as e: