diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 04b7a66..5eb8c78 100644 --- a/lidar_pipeline/dtm.py +++ b/lidar_pipeline/dtm.py @@ -157,7 +157,6 @@ def validate_laz(laz_file): pass # Fallback: try PDAL (handles COPC v1.1 that laspy can't read) - import subprocess try: result = subprocess.run( ["pdal", "info", str(laz_file), "--summary"], diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index ea10369..d4660cc 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -4,7 +4,7 @@ LidarArchaeoPipeline coordinates the full processing chain: 1. Ground classification (PDAL/SMRF) 2. DTM generation 3. Visualization generation (17 products) -4. Rendering (WebP + PDF report) +4. Rendering (AVIF/WebP conversion) """ import logging @@ -61,9 +61,9 @@ from .visualizations import ( generate_lrm, generate_openness, generate_mslrm, generate_tpi, generate_sailore, generate_roughness, generate_wavelet, - generate_svf, generate_aniso_open, + generate_svf, generate_aniso_open, generate_paths, ) -from .gpu import gpu_cleanup, num_gpus +from .gpu import gpu_cleanup, num_gpus, safe_gpu_call from .ign import generate_ign_overlay from .rendering import tif_to_png @@ -85,6 +85,7 @@ VIZ_STEPS = [ ('roughness', generate_roughness), ('svf', generate_svf), ('aniso_open', generate_aniso_open), + ('paths', generate_paths), ('wavelet', generate_wavelet), ('ortho', lambda d, b, v, r: generate_ign_overlay( d, b, v, r, @@ -278,10 +279,11 @@ class LidarArchaeoPipeline: t0 = time.time() try: # IGN overlays don't use SharedDEM (they download external data) + # Non-IGN visualizations use safe_gpu_call for GPU→CPU fallback if name in ('ortho', 'topo'): result = func(dtm_file, basename, file_vis_dir, resolution) else: - result = func(dtm_file, basename, file_vis_dir, resolution, shared=shared) + result = safe_gpu_call(func, dtm_file, basename, file_vis_dir, resolution, shared=shared) vis_results[name] = result elapsed = time.time() - t0 if result: @@ -524,10 +526,6 @@ class LidarArchaeoPipeline: try: if self.temp_dir.exists(): shutil.rmtree(self.temp_dir) - # Also clean up any subdirectories inside temp/ - temp_base = self.output_dir / "temp" - if temp_base.exists(): - shutil.rmtree(temp_base) logger.info(" ✓ Fichiers temporaires supprimés") except Exception as e: logger.warning(f" Note: Impossible de supprimer les fichiers temporaires: {e}") diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index f490345..93ad007 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -80,57 +80,64 @@ _FRANCE_OUTLINE_L93 = np.array([ COLORMAPS = { # === Famille RELIEF : rouge=surélévation, bleu=dépression === - # Roma (Crameri): perceptually uniform, CVD-friendly, dark center → near-zero values visible - # Falls back to RdBu_r if cmcrameri unavailable + # Diverging: rouge vif=positif, bleu vif=négatif, blanc=plat 'curvature': { - 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'cmap': 'bwr', '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), }, 'mslrm': { - 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'cmap': 'seismic', 'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)', 'legend': 'Relief combiné multi-échelles\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nDé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': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'cmap': 'seismic', '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), }, 'tpi': { - 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'cmap': 'seismic', 'title': 'TPI - Topographic Position Index (4 échelles)', 'legend': 'Position dans le paysage\nRouge = Plus haut que le voisinage (crête, plateau)\nBleu = Plus bas que le voisinage (fossé, vallée)\nCombine 4 échelles : 3m, 15m, 50m, 200m', 'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, 'sailore': { - 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'cmap': 'seismic', 'title': 'SAILORE - LRM Auto-Adaptatif', 'legend': 'Relief local adaptatif\nRouge = Surélévation | Bleu = Dépression\n\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), }, 'aniso_open': { - 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', + 'cmap': 'seismic', 'title': 'Openness Anisotropique (pondération directionnelle)', 'legend': 'Openness positive - négative pondérée (degrés)\nRouge = Surélévation dominante (mur, levée)\nBleu = Dépression dominante (fossé, doline)\nPondère les directions NW-SE et NE-SW davantage', 'description': 'Openness avec pondération anisotropique — détecte mieux les structures alignées NW-SE et NE-SW', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, - # === Famille OUVERTURE : clair=ouvert, sombre=fermé === + 'paths': { + 'cmap': 'inferno', + 'title': 'Cheminement (chemins et sentiers)', + 'legend': 'Openness directionnelle maximale\nJaune/vif = Chemin ou sentier\nNoir = Terrain plat\n\nDétecte les structures linéaires dans toutes les directions', + 'description': 'Différence openness positive-négative maximale sur 8 directions — chemins, sentiers, ornières', + 'vmin_mode': 'percentile', 'vmin_pct': 5, + 'vmax_mode': 'percentile', 'vmax_pct': 99, + }, + # === Famille OUVERTURE : séquentiel, toujours positif === '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)', + 'legend': 'Angle d\'ouverture vers le ciel (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, + 'vmin_mode': 'percentile', 'vmin_pct': 2, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, 'negative_openness': { @@ -138,14 +145,14 @@ COLORMAPS = { '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, + 'vmin_mode': 'percentile', 'vmin_pct': 2, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, 'svf': { - 'cmap': 'bone_r', + 'cmap': 'hot_r', 'title': 'Sky-View Factor (fraction de ciel visible)', - 'legend': 'Proportion de ciel visible depuis chaque point\nClair = Ciel dégagé (sommet, plateau, levée)\nSombre = Ciel masqué (vallée, fossé, tranchée)\nContraste adapté aux valeurs réelles (percentiles 2-98)', - 'description': 'Détection de micro-relief — fossés sombres, levées claires, complémentaire de l\'openness', + 'legend': 'Proportion de ciel visible depuis chaque point\nBlanc/jaune = Ciel masqué (vallée, fossé, tranchée)\nNoir = Ciel dégagé (sommet, plateau)\nLes fossés ressortent en vif — excellent pour structures linéaires', + 'description': 'Détection de micro-relief — fossés en jaune/blanc, levées en sombre', 'vmin_mode': 'percentile', 'vmin_pct': 2, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, @@ -155,16 +162,16 @@ COLORMAPS = { '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, + 'vmin_mode': 'percentile', 'vmin_pct': 1, + 'vmax_mode': 'percentile', 'vmax_pct': 99, }, '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', + 'legend': 'Inclinaison en degrés\nMin: {vmin:.1f}° | Max: {vmax:.1f}°\nJaune = Forte pente | Violet foncé = Terrain plat', + 'description': 'Murs, talus et bords ressortent en jaune — terrain plat en sombre', 'vmin_mode': 'fixed', 'vmin_val': 0, - 'vmax_mode': 'percentile', 'vmax_pct': 95, + 'vmax_mode': 'percentile', 'vmax_pct': 97, }, 'aspect': { 'cmap': 'twilight', @@ -175,12 +182,12 @@ COLORMAPS = { 'vmax_mode': 'fixed', 'vmax_val': 360, }, 'roughness': { - 'cmap': 'magma', + 'cmap': 'plasma', 'title': 'Rugosité Multi-Échelle (3m + 15m)', - 'legend': 'Irrégularité du terrain combinée fine + large\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nCombine rugosité fine 3m (70%) + large 15m (30%)', + 'legend': 'Irrégularité du terrain combinée fine + large\nViolet foncé = Surface lisse (route, mur, sol plat)\nJaune vif = Surface rugueuse (végétation, ruines, pierres)\nCombine rugosité fine 3m (70%) + large 15m (30%)', 'description': 'Mesure la variabilité locale — surfaces anthropiques lisses vs naturelles rugueuses', 'vmin_mode': 'fixed', 'vmin_val': 0, - 'vmax_mode': 'percentile', 'vmax_pct': 97, + 'vmax_mode': 'percentile', 'vmax_pct': 98, }, 'wavelet': { 'cmap': 'cividis', @@ -345,12 +352,15 @@ def _nice_scale(extent_m): Returns (scale_m, label) where label is like '100 m' or '500 m' or '1 km'. """ - nice_scales = [50, 100, 200, 500, 1000, 2000, 5000, 10000] - # Pick the largest scale <= 20% of extent + nice_scales = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000] + # Pick the largest scale that fits within 20% of extent + max_scale = extent_m * 0.20 chosen = nice_scales[0] for s in nice_scales: - if s <= extent_m * 0.20: + if s <= max_scale: chosen = s + else: + break if chosen >= 1000: return chosen, f"{chosen // 1000} km" return chosen, f"{chosen} m" @@ -496,10 +506,11 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, im = ax.imshow(data, cmap=cmap, aspect='equal', origin='upper', interpolation='bilinear') - ax.set_title(f"{title}", fontsize=14, fontweight='bold', pad=8) - ax.text(0.5, 1.01, description, transform=ax.transAxes, - fontsize=10, fontstyle='italic', color='#555555', - ha='center', va='bottom') + ax.set_title(f"{title}", fontsize=14, fontweight='bold', pad=10) + if description: + ax.text(0.5, 1.04, description, transform=ax.transAxes, + fontsize=10, fontstyle='italic', color='#555555', + ha='center', va='bottom') # Colorbar/legend area — full height alongside data cbar_left = data_left + data_width_frac + 0.02 @@ -569,34 +580,35 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, spine.set_color('black') spine.set_linewidth(0.8) - # North arrow — compass rose style, inside the data area (top-right corner) + # North arrow — compass rose in bottom-right corner of data area # Semi-transparent background for readability over any data - north_ax = fig.add_axes([data_left + data_width_frac - 0.06, - data_bottom + data_height_frac - 0.12, - 0.05, 0.10], + north_ax = fig.add_axes([data_left + data_width_frac - 0.07, + data_bottom + 0.01, + 0.06, 0.14], facecolor='none') - north_ax.set_xlim(-1.2, 1.2) - north_ax.set_ylim(-0.3, 1.5) + north_ax.set_xlim(-1.5, 1.5) + north_ax.set_ylim(-1.5, 1.5) north_ax.axis('off') north_ax.set_aspect('equal') + # Compass rose centered at (0, 0) — all 4 cardinals equidistant from center # Semi-transparent white background circle - circle_bg = plt.Circle((0, 0.5), 0.85, facecolor='white', edgecolor='#888888', + circle_bg = plt.Circle((0, 0), 1.0, facecolor='white', edgecolor='#888888', linewidth=0.5, alpha=0.7, zorder=1) north_ax.add_patch(circle_bg) - # N arrow - north_ax.annotate('N', xy=(0, 1.3), fontsize=9, fontweight='bold', + # N arrow (pointing up = North) + north_ax.annotate('N', xy=(0, 1.35), fontsize=9, fontweight='bold', ha='center', va='bottom', color='#b22222', zorder=10) - north_ax.plot([0, 0], [0.0, 1.0], color='#b22222', linewidth=2.0, zorder=10) - north_ax.add_patch(MplPolygon([[0, 0.3], [-0.2, 0.7], [0, 1.0], [0.2, 0.7]], + north_ax.plot([0, 0], [-0.5, 1.0], color='#b22222', linewidth=2.0, zorder=10) + north_ax.add_patch(MplPolygon([[0, 0.5], [-0.2, 0.7], [0, 1.0], [0.2, 0.7]], closed=True, facecolor='#b22222', edgecolor='#b22222', zorder=9)) - # Cardinal ticks - for angle, label in [(90, ''), (0, 'E'), (180, 'O'), (270, 'S')]: + # Cardinal ticks — all centered at (0, 0) + for angle, label in [(90, 'N'), (0, 'E'), (180, 'O'), (270, 'S')]: rad = np.radians(angle) - north_ax.plot([0.85*np.cos(rad), 1.05*np.cos(rad)], - [0.85*np.sin(rad), 1.05*np.sin(rad)], + north_ax.plot([1.0*np.cos(rad), 1.2*np.cos(rad)], + [1.0*np.sin(rad), 1.2*np.sin(rad)], color='#555555', linewidth=0.8, zorder=5) if label: - north_ax.text(1.15*np.cos(rad), 1.15*np.sin(rad), label, + north_ax.text(1.35*np.cos(rad), 1.35*np.sin(rad), label, fontsize=6, ha='center', va='center', color='#555555', zorder=5) # Bottom info bar — enriched with source, method, date @@ -621,7 +633,9 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, else: line1_parts.append(f"X: {min_x:.0f}–{max_x:.0f} Y: {min_y:.0f}–{max_y:.0f}") line1_parts.append(f"EPSG:2154") - line1_parts.append(f"Res: {resolution}m/px") + # Round resolution to avoid ugly decimals like 0.499999959 + res_display = round(resolution, 2) if resolution < 1 else round(resolution, 1) + line1_parts.append(f"Res: {res_display}m/px") line1_parts.append(f"Emprise: {extent_km_x:.1f}×{extent_km_y:.1f}km") if not is_rgb: line1_parts.append(f"Alt: {alt_min:.1f}–{alt_max:.1f}m") @@ -793,15 +807,24 @@ def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): 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 + # Look for images in per-file subdirectory first, then fallback to main dir file_vis_dir = vis_dir / basename + png_files = [] if file_vis_dir.exists(): - png_files = sorted(file_vis_dir.glob("*.webp")) + png_files = sorted(file_vis_dir.glob("*.avif")) + sorted(file_vis_dir.glob("*.webp")) else: - png_files = sorted(vis_dir.glob(f"{basename}_*.webp")) - if not png_files: + png_files = sorted(vis_dir.glob(f"{basename}_*.avif")) + sorted(vis_dir.glob(f"{basename}_*.webp")) + # Deduplicate in case both formats exist + seen = set() + unique_files = [] + for f in png_files: + if f not in seen: + seen.add(f) + unique_files.append(f) + if not unique_files: logger.warning(f" ✗ Aucune image trouvée pour {basename}") return None + png_files = unique_files # Categorize situ_files = [] diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index f325736..96acd10 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -15,19 +15,41 @@ from pathlib import Path import numpy as np import rasterio -from scipy.ndimage import generic_filter -from scipy.stats import binned_statistic_2d - from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, xp_minimum_filter, xp_maximum_filter, gpu_cleanup +from . import gpu as _gpu_mod logger = logging.getLogger("lidar") -# Use CuPy array module when available -if HAS_GPU: - import cupy as cp - xp = cp -else: - xp = np +# CuPy module reference — lazily imported on first GPU use. +# If disable_gpu() is called at runtime, HAS_GPU becomes False +# and xp delegates to numpy instead. +_cp = None + + +class _XPProxy: + """Proxy that delegates array operations to cupy or numpy. + + Checks HAS_GPU on every attribute access so that disable_gpu() + (called on CUDA errors) takes effect immediately, without needing + to change every call site in visualizations.py. + """ + + def __getattr__(self, name): + global _cp + from . import gpu as _gpu_mod + if _gpu_mod.HAS_GPU: + if _cp is None: + try: + import cupy + _cp = cupy + except ImportError: + pass + if _cp is not None: + return getattr(_cp, name) + return getattr(np, name) + + +xp = _XPProxy() class SharedDEM: @@ -118,14 +140,14 @@ class SharedDEM: @property def filled_gpu(self): """Lazy GPU copy of the filled DEM.""" - if self._filled_gpu is None and HAS_GPU: + if self._filled_gpu is None and _gpu_mod.HAS_GPU: self._filled_gpu = to_gpu(self.filled) return self._filled_gpu @property def dem_gpu(self): """Lazy GPU copy of the DEM.""" - if self._dem_gpu is None and HAS_GPU: + if self._dem_gpu is None and _gpu_mod.HAS_GPU: self._dem_gpu = to_gpu(self.dem_np) return self._dem_gpu @@ -134,10 +156,14 @@ def _filter_nanaware_from_filled(shared, filter_func, *args, **kwargs): """Apply filter on pre-filled DEM data (skips expensive _fill_nans). Uses the SharedDEM.filled array directly, then restores NaN mask. - If GPU is available, uses the lazy GPU copy to avoid CPU↔GPU transfers. + If GPU is available, reuses the lazy GPU copy to avoid redundant transfers. """ - if HAS_GPU and shared.filled_gpu is not None: - filled_gpu = to_gpu(shared.filled) + if _gpu_mod.HAS_GPU: + filled_gpu = shared.filled_gpu + else: + filled_gpu = None + + if filled_gpu is not None: result_gpu = filter_func(filled_gpu, *args, **kwargs) result = to_cpu(result_gpu) gpu_cleanup() @@ -227,15 +253,16 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs): Returns: Filtered array with original NaN positions preserved. """ - is_gpu_arr = HAS_GPU and isinstance(arr, cp.ndarray) + is_gpu_arr = _gpu_mod.HAS_GPU and _cp is not None and isinstance(arr, _cp.ndarray) arr_np = to_cpu(arr) if is_gpu_arr else arr filled, nan_mask = _fill_nans(arr_np) - if use_gpu and HAS_GPU: + if use_gpu and _gpu_mod.HAS_GPU: filled_gpu = to_gpu(filled) result_gpu = filter_func(filled_gpu, *args, **kwargs) result = to_cpu(result_gpu) + gpu_cleanup() else: result = filter_func(filled, *args, **kwargs) @@ -254,7 +281,7 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): Applies percentile normalization and gamma correction to restore contrast lost by averaging multiple azimuths. """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Hillshade multidirectionnel{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_hillshade_multi.tif" @@ -264,9 +291,9 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): transform = shared.transform crs = shared.crs dem = to_gpu(shared.dem_np) - dy = to_gpu(shared.dy) if HAS_GPU else shared.dy - dx = to_gpu(shared.dx) if HAS_GPU else shared.dx - slope = to_gpu(shared.slope_rad) if HAS_GPU else shared.slope_rad + dy = to_gpu(shared.dy) if _gpu_mod.HAS_GPU else shared.dy + dx = to_gpu(shared.dx) if _gpu_mod.HAS_GPU else shared.dx + slope = to_gpu(shared.slope_rad) if _gpu_mod.HAS_GPU else shared.slope_rad aspect = xp.arctan2(dy, dx) sin_slope = xp.sin(slope) cos_slope = xp.cos(slope) @@ -299,7 +326,7 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): # Contrast enhancement: percentile stretch + gamma combined_np = to_cpu(combined) - nan_mask = shared.nan_mask if shared else np.isnan(to_cpu(dem_np) if HAS_GPU else dem_np) + nan_mask = shared.nan_mask if shared else np.isnan(dem_np) valid = combined_np[~nan_mask] if len(valid) > 0: p2, p98 = np.percentile(valid, 2), np.percentile(valid, 98) @@ -319,7 +346,7 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): """Generate slope map (degrees) — GPU if available.""" - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Pente (Slope){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_slope.tif" @@ -330,7 +357,7 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): crs = shared.crs slope = shared.slope_deg nan_mask = shared.nan_mask - if HAS_GPU: + if _gpu_mod.HAS_GPU: slope = to_gpu(slope) else: dem_np, transform, crs = _read_dem(dem_file) @@ -338,7 +365,7 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): dy, dx = xp.gradient(dem) slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) * 180 / xp.pi nan_mask = np.isnan(dem_np) - _save_tif(output, to_cpu(slope) if HAS_GPU else slope, transform, crs, nan_mask=nan_mask) + _save_tif(output, to_cpu(slope) if _gpu_mod.HAS_GPU else slope, transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -348,7 +375,7 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): """Generate aspect (slope orientation) map — GPU if available.""" - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Aspect (Orientation){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_aspect.tif" @@ -359,7 +386,7 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): crs = shared.crs aspect = shared.aspect nan_mask = shared.nan_mask - if HAS_GPU: + if _gpu_mod.HAS_GPU: aspect = to_gpu(aspect) else: dem_np, transform, crs = _read_dem(dem_file) @@ -368,7 +395,7 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): aspect = xp.arctan2(dy, dx) * 180 / xp.pi aspect = xp.mod(aspect, 360) nan_mask = np.isnan(dem_np) - _save_tif(output, to_cpu(aspect) if HAS_GPU else aspect, transform, crs, nan_mask=nan_mask) + _save_tif(output, to_cpu(aspect) if _gpu_mod.HAS_GPU else aspect, transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -378,7 +405,7 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): """Generate curvature (terrain concavity/convexity) map — GPU if available.""" - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Courbure (Curvature){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_curvature.tif" @@ -390,7 +417,7 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): dx = shared.dx dy = shared.dy nan_mask = shared.nan_mask - if HAS_GPU: + if _gpu_mod.HAS_GPU: dx = to_gpu(dx) dy = to_gpu(dy) else: @@ -420,7 +447,7 @@ def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None): Kernel sigma adapts to resolution: finer kernel at higher resolution to capture micro-relief details. At 0.5m/px: 15m, at 0.2m/px: ~5m. """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Local Relief Model{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_lrm.tif" @@ -454,7 +481,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): 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 "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Sky-View Factor (ray-tracing){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_svf.tif" @@ -466,7 +493,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution - dem = to_gpu(shared.filled) if HAS_GPU else shared.filled + dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) @@ -474,7 +501,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if HAS_GPU else filled + dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 16 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) @@ -509,7 +536,8 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): 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 uses cos²(horizon angle) — fraction of visible sky + svf += xp.cos(horizon) ** 2 svf /= n_dirs svf_np = to_cpu(svf).astype(np.float32) @@ -532,7 +560,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh Ray radius adapts to resolution: 100m for better detection of large enclosures. """ name = "positive_openness" if positive else "negative_openness" - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_{name}.tif" @@ -544,7 +572,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution - dem = to_gpu(shared.filled) if HAS_GPU else shared.filled + dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) @@ -552,7 +580,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if HAS_GPU else filled + dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) @@ -597,71 +625,13 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh return None -def generate_local_dominance(dem_file, basename, vis_dir, resolution, shared=None, - radius=15, pmin=2, pmax=98): - """Local Dominance — proportion of neighborhood below center point. - - LD = (dem - local_min) / (local_max - local_min + epsilon) - - High values = locally dominant (peak, ridge) - Low values = locally recessed (valley, pit) - - Uses minimum/maximum filters on the filled DEM, then restores NaN mask. - Complements openness by measuring local height position rather than angular extent. - """ - gpu_tag = " [GPU]" if HAS_GPU else "" - logger.info(f" → Dominance Locale (rayon {radius}m){gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_local_dominance.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - nan_mask = shared.nan_mask - dem_np = shared.dem_np - else: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) - - radius_px = max(1, int(radius / resolution)) - if radius_px % 2 == 0: - radius_px += 1 - - local_min = _filter_nanaware_from_filled( - shared, xp_minimum_filter, size=radius_px - ) if shared else _filter_nanaware( - dem_np, xp_minimum_filter, size=radius_px - ) - - local_max_data = _filter_nanaware_from_filled( - shared, xp_maximum_filter, size=radius_px - ) if shared else _filter_nanaware( - dem_np, xp_maximum_filter, size=radius_px - ) - - # Local dominance ratio - epsilon = 0.01 # Avoid division by zero on flat terrain - local_range = local_max_data - local_min + epsilon - dominance = (dem_np - local_min) / local_range - dominance = np.clip(dominance, 0, 1) - dominance[nan_mask] = np.nan - - _save_tif(output, dominance.astype(np.float32), transform, crs, nan_mask=nan_mask) - logger.info(f" ✓ Dominance Locale terminée ({time.time()-t0:.1f}s){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur local_dominance: {e}", exc_info=True) - return None - - def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None): """Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available). Scales adapt to resolution. Std normalization per scale. Weighted combination favoring archaeologically relevant scales (5-25m). """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_mslrm.tif" @@ -731,7 +701,7 @@ def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None): Computed at 4 scales with std normalization and weighted combination. Weights favor fine and medium scales (archaeologically relevant). """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → TPI multi-échelle{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_tpi.tif" @@ -796,7 +766,7 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None): Kernel size adapts to local slope: flat areas get larger kernels, steep areas get smaller kernels. Scales adapt to resolution. """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → SAILORE (LRM adaptatif){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_sailore.tif" @@ -816,10 +786,9 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None): slope_deg = np.degrees(slope) slope_deg[nan_mask] = np.nan - # Adaptive scales: finer at higher resolution - sigma_min_m = max(1.0, 2.0 * 0.5 / resolution) # 2m at 0.5, ~5m at 0.2 - sigma_mid_m = max(5.0, 13.5 * 0.5 / resolution) # 13.5m at 0.5, ~33m at 0.2 - sigma_max_m = max(5.0, 25.0 * 0.5 / resolution) # 25m at 0.5, ~62m at 0.2 + # Fixed physical scales (independent of resolution) + sigma_min_m = 2.0 # 2m — fine detail + sigma_max_m = 25.0 # 25m — broad relief sigma_min = sigma_min_m / resolution sigma_max = sigma_max_m / resolution sigma_mid = (sigma_min + sigma_max) / 2 @@ -870,7 +839,7 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): Combines fine (3m) and broad (15m) roughness for better detection of archaeological features at multiple scales. """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Rugosité de surface{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_roughness.tif" @@ -924,7 +893,6 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): roughness = 0.7 * roughness_fine / fine_std + 0.3 * roughness_broad / broad_std roughness[nan_mask] = np.nan - roughness = to_cpu(roughness) _save_tif(output, roughness, transform, crs) logger.info(f" ✓ Rugosité terminée ({time.time()-t0:.1f}s){gpu_tag}") return output @@ -933,90 +901,6 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): return None -# ============================================================ -# Anomalies -# ============================================================ - -def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None): - """Statistical anomaly detection - std-normalized multi-scale relief + Local Moran's I — GPU if available. - - Uses MSRM (multi-scale LRM) instead of single-scale LRM for better detection - of anomalies at all scales. - """ - gpu_tag = " [GPU]" if HAS_GPU else "" - logger.info(f" → Détection anomalies statistiques{gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_anomalies.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - dem_np = shared.dem_np - nan_mask = shared.nan_mask - else: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) - - # Multi-scale LRM: compute MSRM-like combined relief - min_scale = max(2.0, resolution * 4) - candidate_scales = [2, 5, 10, 20, 50, 100] - sigmas = [s for s in candidate_scales if s >= min_scale] - lrm_stack = [] - - for sigma in sigmas: - sigma_px = sigma / resolution - if shared: - local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px) - else: - local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px) - lrm = dem_np - local_mean - lrm[nan_mask] = np.nan - # Std normalization — preserves contrast better than z-score - valid_lrm = lrm[~nan_mask] - lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 - lrm_norm = lrm / lrm_std - lrm_stack.append(lrm_norm.astype(np.float32)) - - # Weighted RMS combination (favor 5-25m scales) - scale_weights = {2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6} - weights = np.array([scale_weights.get(s, 1.0) for s in sigmas]) - lrm_array = np.array(lrm_stack) - weights_3d = weights[:, np.newaxis, np.newaxis] - with np.errstate(invalid='ignore', divide='ignore'): - with warnings.catch_warnings(): - warnings.filterwarnings('ignore', message='Mean of empty slice') - msrm = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights)) - msrm[nan_mask] = np.nan - - # Std normalization of MSRM — preserves contrast better than z-score - valid_msrm = msrm[~nan_mask] - msrm_std = max(np.nanstd(valid_msrm), 0.01) if len(valid_msrm) > 0 else 0.01 - z_score = msrm / msrm_std - - # Local Moran's I for spatial clustering - window = max(3, int(10 / resolution)) - if window % 2 == 0: - window += 1 - - if shared: - local_mean_z = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window) - else: - local_mean_z = _filter_nanaware(z_score, xp_uniform_filter, size=window) - z_mean_global = np.nanmean(z_score[~nan_mask]) if np.any(~nan_mask) else 0 - z_std_global = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01 - morans_i = z_score * (local_mean_z - z_mean_global) / z_std_global - anomaly_score = np.abs(z_score) * np.sign(morans_i) - anomaly_score[nan_mask] = np.nan - - _save_tif(output, anomaly_score.astype(np.float32), transform, crs) - logger.info(f" ✓ Anomalies terminé ({time.time()-t0:.1f}s){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur anomalies: {e}", exc_info=True) - return None - - # ============================================================ # Wavelet # ============================================================ @@ -1032,7 +916,7 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): Uses std normalization per scale and weighted combination with emphasis on archaeologically relevant scales (2-50m). """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Ondelette Mexican Hat multi-échelle{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_wavelet.tif" @@ -1072,10 +956,14 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): 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(to_gpu(filled), sigma=sigma_px) - response = to_cpu(response) + if _gpu_mod.HAS_GPU: + try: + from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace + response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px) + response = to_cpu(response) + except Exception: + from scipy.ndimage import gaussian_laplace + response = -gaussian_laplace(filled, sigma=sigma_px) else: from scipy.ndimage import gaussian_laplace response = -gaussian_laplace(filled, sigma=sigma_px) @@ -1106,200 +994,28 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): # ============================================================ -# Flow accumulation +# Anisotropic Openness +# ============================================================ +# Path Detection (chemins et sentiers) # ============================================================ -def _d8_accumulate_numba(flow_dir, nodata_mask, rows, cols): - """JIT-compiled D8 flow accumulation loop. +def generate_paths(dem_file, basename, vis_dir, resolution, shared=None): + """Cheminement — openness directionnelle maximale pour détecter chemins et sentiers. - Uses numba for ~100x speedup over pure Python loop. - Falls back to pure Python if numba is unavailable. + Pour chaque direction (8 directions), calcule openness positive - négative, + puis prend le maximum sur toutes les directions. Les chemins et sentiers + ressortent en valeurs élevées quelle que soit leur orientation. + + Contrairement à l'openness anisotropique qui privilégie NW-SE et NE-SW, + cette visualisation traite toutes les directions de manière égale et + combine positive et négative en une seule image. Les chemins perpendiculaires + à une direction auront une forte différence dans cette direction, donc + le maximum sur toutes les directions les fait ressortir. """ - try: - from numba import njit - - @njit(cache=True) - def _accumulate(flow_dir, nodata_mask, rows, cols): - dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int8) - dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int8) - - flow_acc = np.ones((rows, cols), dtype=np.float32) - - # Sort cells by elevation (high to low) — walk downhill - # We use the fact that flow_dir already encodes steepest descent - # Process from highest to lowest elevation - for r in range(rows): - for c in range(cols): - if nodata_mask[r, c]: - flow_acc[r, c] = 0.0 - continue - - # Iterative accumulation: process cells in top-down order - # Multiple passes until convergence - for _pass in range(10): - changed = 0 - for r in range(rows): - for c in range(cols): - if nodata_mask[r, c]: - continue - d = flow_dir[r, c] - if d < 0: - continue - nr = r + dy8[d] - nc = c + dx8[d] - if 0 <= nr < rows and 0 <= nc < cols and not nodata_mask[nr, nc]: - old_acc = flow_acc[nr, nc] - flow_acc[nr, nc] += flow_acc[r, c] - if flow_acc[nr, nc] != old_acc: - changed += 1 - if changed == 0: - break - - return flow_acc - - return _accumulate(flow_dir, nodata_mask, rows, cols) - - except ImportError: - # Fallback: pure Python - return None - - -def _priority_flood(dem, nodata_mask): - """Priority-flood algorithm for sink filling (Wang & Liu 2006). - - O(n log n) compared to 50 iterations of minimum_filter. - Fills pits so water can flow downhill. - """ - import heapq - - rows, cols = dem.shape - filled = dem.copy() - closed = nodata_mask.copy() - open_queue = [] - - # Initialize border cells - for r in range(rows): - for c in [0, cols - 1]: - if not closed[r, c]: - heapq.heappush(open_queue, (filled[r, c], r, c)) - closed[r, c] = True - for c in range(1, cols - 1): - for r in [0, rows - 1]: - if not closed[r, c]: - heapq.heappush(open_queue, (filled[r, c], r, c)) - closed[r, c] = True - - dx8 = [1, 1, 0, -1, -1, -1, 0, 1] - dy8 = [0, 1, 1, 1, 0, -1, -1, -1] - - while open_queue: - elev, r, c = heapq.heappop(open_queue) - for d in range(8): - nr, nc = r + dy8[d], c + dx8[d] - if 0 <= nr < rows and 0 <= nc < cols and not closed[nr, nc]: - if filled[nr, nc] < elev: - filled[nr, nc] = elev # Fill the pit - closed[nr, nc] = True - heapq.heappush(open_queue, (filled[nr, nc], nr, nc)) - - return filled - - -def generate_flow(dem_file, basename, vis_dir, resolution, shared=None): - """Flow accumulation using D8 algorithm — priority-flood sink filling, accumulation via numba.""" - gpu_tag = " [GPU]" if HAS_GPU else "" - logger.info(f" → Accumulation de flux D8{gpu_tag}...") + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" + logger.info(f" → Cheminement (chemins et sentiers){gpu_tag}...") t0 = time.time() - output = vis_dir / f"{basename}_flow.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - dem_np = shared.dem_np - nodata_mask = shared.nan_mask - else: - dem_np, transform, crs = _read_dem(dem_file) - nodata_mask = np.isnan(dem_np) - - rows, cols = dem_np.shape - - # Sink filling — priority-flood (O(n log n), faster than 50× minimum_filter) - dem_filled_np = _priority_flood(dem_np, nodata_mask) - - # D8 slope — vectorized - dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int32) - dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int32) - dist8 = np.array([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.zeros((rows, cols), dtype=np.float64) - - padded = np.pad(dem_filled_np, 1, mode='constant', - constant_values=np.nanmax(dem_filled_np[~np.isnan(dem_filled_np)]) + 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_np - neighbor_elev) / (dist8[d] * resolution) - slope[nodata_mask] = -1 - better = slope > max_slope - flow_dir[better] = d - max_slope[better] = slope[better] - - # D8 accumulation — try numba first, fallback to Python - result = _d8_accumulate_numba(flow_dir, nodata_mask.astype(np.bool_), rows, cols) - - if result is not None: - flow_acc = result - logger.info(f" Accumulation D8 via numba") - else: - # Pure Python fallback (slow for large DEMs) - logger.info(f" Accumulation D8 via Python (installez numba pour accélérer)") - flat_dem = dem_filled_np[~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){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur flux: {e}", exc_info=True) - return None - - -# ============================================================ -# Sky-View Factor (SVF) -# ============================================================ - -def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): - """Sky-View Factor - fraction of sky visible from each point (GPU if available). - - SVF = average of cos²(horizon_angle) across 16 directions. - High SVF (near 1) = open sky (ridgetop, plateau) - Low SVF (near 0) = enclosed sky (valley, deep trench) - - Excellent for detecting archaeological earthworks: ditches appear dark, - embankments appear bright. Complements openness which uses raw angles. - """ - gpu_tag = " [GPU]" if HAS_GPU else "" - logger.info(f" → Sky-View Factor{gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_svf.tif" + output = vis_dir / f"{basename}_paths.tif" try: if shared: @@ -1308,7 +1024,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution - dem = to_gpu(shared.filled) if HAS_GPU else shared.filled + dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) @@ -1316,21 +1032,24 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if HAS_GPU else filled + dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled - n_dirs = 16 # More directions for smoother SVF + n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx_dir = np.cos(angles) dy_dir = np.sin(angles) max_dist = min(int(100 / res), 300) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - svf_sum = xp.zeros_like(dem) + max_diff = xp.full_like(dem, -1e6) # Will track max over all directions for d_idx in range(n_dirs): ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] - # Find maximum horizon elevation angle in this direction - max_horizon_angle = xp.zeros_like(dem) + + # Positive openness: max zenith angle in this direction + max_pos_angle = xp.zeros_like(dem) + # Negative openness: max nadir angle in this direction + max_neg_angle = xp.zeros_like(dem) for step in range(1, max_dist + 1): px = int(round(ddx * step)) @@ -1342,27 +1061,30 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): elev_diff = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] - dem - # Horizon angle from horizontal (positive = terrain above viewer) - angle = xp.arctan2(elev_diff, dist_m) - max_horizon_angle = xp.where(xp.isnan(angle), max_horizon_angle, - xp.maximum(max_horizon_angle, xp.nan_to_num(angle, nan=0))) + # Positive: angle to terrain above viewer + pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) + max_pos_angle = xp.where(xp.isnan(pos_angle), max_pos_angle, + xp.maximum(max_pos_angle, xp.nan_to_num(pos_angle, nan=0))) - # SVF uses cos²(horizon angle) — fraction of visible sky in this direction - cos2 = xp.cos(max_horizon_angle) ** 2 - svf_sum += cos2 + # Negative: angle to terrain below viewer + neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) + max_neg_angle = xp.where(xp.isnan(neg_angle), max_neg_angle, + xp.maximum(max_neg_angle, xp.nan_to_num(neg_angle, nan=0))) - svf_result = to_cpu(svf_sum / n_dirs).astype(np.float32) - svf_result[nan_mask] = np.nan - _save_tif(output, svf_result, transform, crs) - logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}") + # Difference highlights linear features perpendicular to this direction + diff = max_pos_angle - max_neg_angle + max_diff = xp.maximum(max_diff, diff) + + paths_result = to_cpu(max_diff).astype(np.float32) + paths_result[nan_mask] = np.nan + _save_tif(output, paths_result, transform, crs) + logger.info(f" ✓ Cheminement terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: - logger.error(f" ✗ Erreur SVF: {e}", exc_info=True) + logger.error(f" ✗ Erreur cheminement: {e}", exc_info=True) return None -# ============================================================ -# Anisotropic Openness # ============================================================ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None): @@ -1375,7 +1097,7 @@ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None): The anisotropic weighting makes subtle linear features more visible than standard isotropic openness which averages all directions equally. """ - gpu_tag = " [GPU]" if HAS_GPU else "" + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Openness Anisotropique{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_aniso_open.tif" @@ -1387,7 +1109,7 @@ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None): dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution - dem = to_gpu(shared.filled) if HAS_GPU else shared.filled + dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) @@ -1395,7 +1117,7 @@ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None): res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if HAS_GPU else filled + dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)