"""Terrain visualization functions for LiDAR archaeological analysis. Each function takes (dem_file, basename, vis_dir, resolution) as explicit parameters and returns the path to the output GeoTIFF file, or None on error. When a SharedDEM object is provided via the `shared` parameter, pre-computed data (gradient, NaN mask, LRM) is reused across visualizations to avoid redundant I/O and computation. """ import logging import math import time import warnings from pathlib import Path import numpy as np import rasterio from .gpu import to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, gpu_cleanup from . import gpu as _gpu_mod logger = logging.getLogger("lidar") # 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 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: """Pre-computed DEM data shared across all visualizations. Reads the DEM once and lazily computes on first access: - NaN mask and filled DEM (avoids 20+ calls to _fill_nans) - Gradient components (shared by hillshade, slope) - LRM at 15m kernel (shared by mslrm + sailore) Attributes are computed lazily on first access to avoid computing data that is never used (e.g. LRM when only hillshade needs generation). """ def __init__(self, dem_file, resolution): dem_np, transform, crs = _read_dem(dem_file) self.dem_file = dem_file self.resolution = resolution self.transform = transform self.crs = crs self.nan_mask = np.isnan(dem_np) self.dem_np = dem_np.astype(np.float32) # Lazy caches — computed on first access self._filled = None self._gradient = None # (dy, dx, slope_rad, slope_deg) self._lrm_15 = None # GPU lazy caches self._filled_gpu = None self._dem_gpu = None @property def filled(self): """Filled DEM (NaN interpolated) — computed lazily.""" if self._filled is None: logger.debug(" → Calcul filled DEM (interpolation NaN)...") self._filled, _ = _fill_nans(self.dem_np) return self._filled @property def dy(self): self._ensure_gradient() return self._gradient[0] @property def dx(self): self._ensure_gradient() return self._gradient[1] @property def slope_rad(self): self._ensure_gradient() return self._gradient[2] @property def slope_deg(self): self._ensure_gradient() return self._gradient[3] @property def lrm_15(self): """LRM at 15m kernel — computed lazily.""" if self._lrm_15 is None: logger.debug(" → Calcul LRM 15m...") sigma_15 = 15.0 / self.resolution local_mean_15 = _filter_nanaware_from_filled(self, xp_gaussian_filter, sigma=sigma_15) self._lrm_15 = self.dem_np - local_mean_15 self._lrm_15[self.nan_mask] = np.nan return self._lrm_15 def _ensure_gradient(self): """Compute gradient components lazily on first access.""" if self._gradient is None: logger.debug(" → Calcul gradient...") dy = dx = None if _gpu_mod.is_gpu_active() and self.filled_gpu is not None: try: g = self.filled_gpu dy = to_cpu(xp.gradient(g, self.resolution, axis=0)) dx = to_cpu(xp.gradient(g, self.resolution, axis=1)) except Exception as e: logger.warning(f" Gradient GPU impossible ({e}) — repli CPU") dy = dx = None if dy is None: dy = np.gradient(self.filled, self.resolution, axis=0) dx = np.gradient(self.filled, self.resolution, axis=1) slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) slope_deg = np.degrees(slope_rad) self._gradient = (dy, dx, slope_rad, slope_deg) @property def filled_gpu(self): """Lazy GPU copy of the filled DEM.""" 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 _gpu_mod.HAS_GPU: self._dem_gpu = to_gpu(self.dem_np) return self._dem_gpu 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, reuses the lazy GPU copy to avoid redundant transfers. """ 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() else: result = filter_func(shared.filled, *args, **kwargs) result[shared.nan_mask] = np.nan return result def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None, nan_mask=None): """Helper to save a 2D or 3D array as GeoTIFF. Args: nan_mask: Optional boolean mask (True=NaN) to apply before saving. Restores NaN zones in gradient-derived products that were computed on the filled DEM. """ if nan_mask is not None: data = np.array(data, dtype=dtype, copy=True) data[nan_mask] = np.nan # Auto-detect nodata for float types with NaN if nodata is None and dtype.startswith('float') and np.any(np.isnan(data)): nodata = float('nan') if data.ndim == 2: height, width = data.shape with rasterio.open( output_path, 'w', driver='GTiff', height=height, width=width, count=count, dtype=dtype, crs=crs, transform=transform, compress='deflate', nodata=nodata ) as dst: dst.write(data.astype(dtype), 1) elif data.ndim == 3: bands, height, width = data.shape with rasterio.open( output_path, 'w', driver='GTiff', height=height, width=width, count=bands, dtype=dtype, crs=crs, transform=transform, compress='deflate', nodata=nodata ) as dst: for i in range(bands): dst.write(data[i].astype(dtype), i + 1) def _read_dem(dem_file): """Read DEM file and return (data, transform, crs).""" with rasterio.open(dem_file) as src: return src.read(1), src.transform, src.crs def _fill_nans(arr): """Fill NaN values using nearest-neighbor interpolation. Returns (filled_array, nan_mask) so the caller can restore NaN after filtering. Via transformée de distance (O(n), vectorisé) : les indices du plus proche voisin valides sortent en une passe. NearestNDInterpolator construisait un cKDTree sur TOUS les points valides (25 M à 0,2 m) — plusieurs secondes par dalle trouée, payées au premier accès de SharedDEM.filled. """ nan_mask = np.isnan(arr) if not np.any(nan_mask): return arr, nan_mask # GPU (cupyx) : la transformée de distance sur 25 M de pixels coûte ~3 s # sur CPU, l'étape la plus longue de la préparation une fois le reste sur GPU. if _gpu_mod.is_gpu_active(): try: from cupyx.scipy.ndimage import distance_transform_edt as edt_gpu cp = _gpu_mod._cp _, idx = edt_gpu(cp.asarray(nan_mask), return_distances=False, return_indices=True) iy, ix = cp.asnumpy(idx[0]), cp.asnumpy(idx[1]) del idx gpu_cleanup() filled = arr.copy() filled[nan_mask] = arr[iy[nan_mask], ix[nan_mask]] return filled, nan_mask except Exception as e: logger.warning(f" Comblement GPU impossible ({e}) — repli CPU") from scipy.ndimage import distance_transform_edt _, (iy, ix) = distance_transform_edt(nan_mask, return_indices=True) filled = arr.copy() filled[nan_mask] = arr[iy[nan_mask], ix[nan_mask]] return filled, nan_mask def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs): """Apply a filter to an array while preserving NaN zones. 1. Fill NaN with nearest-neighbor interpolation 2. Apply the filter 3. Restore original NaN mask on the result Args: arr: Input array (numpy or cupy). filter_func: Function that takes (array, *args, **kwargs) and returns filtered array. use_gpu: If True, apply filter on GPU (send filled array to GPU first). Returns: Filtered array with original NaN positions preserved. """ 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 _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) result[nan_mask] = np.nan return result # ============================================================ # Shared ray-tracing core # ============================================================ def _prepare_dem_for_raycast(dem_file, shared, resolution): """Load DEM and prepare padded array for ray-tracing. Returns (dem_filled, dem_np, rows, cols, res, nan_mask, transform, crs) ready for ray-tracing. dem_filled is a CPU numpy array (filled, no NaN). """ if shared: dem_np = shared.dem_np nan_mask = shared.nan_mask transform = shared.transform crs = shared.crs dem = shared.filled else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) dem = filled res = resolution rows, cols = dem_np.shape return dem, dem_np, rows, cols, res, nan_mask, transform, crs def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=None): """Core ray-tracing: compute max zenith/nadir angles per direction and radius. For each pixel, in each direction, traces rays outward up to max_dist steps, recording the max upward angle (positive openness) and max downward angle (negative openness) reached at each radius checkpoint. Optimisation : on accumule la TANGENTE de l'angle (dz/dist) au lieu de l'angle lui-même — atan étant strictement croissante, max(angles) = atan(max(tangentes)). L'arctan (coûteuse, pleine image) n'est donc plus appliquée qu'aux checkpoints de rayon, pas à chaque pas de rayon. Un seul couple de max cumulés est maintenu, snapshoté à chaque checkpoint (les rayons étant emboîtés, chaque checkpoint réutilisait avant le même calcul 3 fois). fmax ignore les NaN du padding : plus de nan_to_num/where. Padding on CPU (numpy) to avoid GPU memory pressure and pre-compiled kernel mismatches (CUDA_ERROR_NO_BINARY_FOR_GPU on sm_89). The padded array is transferred to GPU once, then each direction is processed and results are streamed back to CPU. Args: dem: CPU numpy array — filled DEM (no NaN), shape (rows, cols). rows, cols: dimensions. res: resolution in m/px. n_dirs: number of directions. max_dist: max ray steps. radii_m: list of radii in meters to record checkpoints. If None, records only at max_dist. Returns: pos_angles: array of shape (n_dirs, n_radii, rows, cols) — max zenith angles neg_angles: array of shape (n_dirs, n_radii, rows, cols) — max nadir angles """ angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx_dir = np.cos(angles) dy_dir = np.sin(angles) if radii_m is not None: radii_steps = [min(int(r / res), max_dist) for r in radii_m] n_radii = len(radii_m) else: radii_steps = [max_dist] n_radii = 1 # Pad on CPU (numpy) — avoids GPU memory pressure and # pre-compiled kernel issues (NO_BINARY_FOR_GPU on sm_89). padded_np = np.pad(dem, max_dist, mode='constant', constant_values=np.nan) # Transfer padded DEM to GPU for computation padded = to_gpu(padded_np) # GPU view of central region — reference elevation for ray-tracing dem = padded[max_dist:max_dist+rows, max_dist:max_dist+cols] # Free the CPU copy — we don't need it anymore del padded_np # Checkpoints triés par pas : (step, r_idx). Les snapshots sont pris quand # le pas courant atteint le pas du checkpoint — les rayons ne dépassent # donc pas le plus grand checkpoint demandé (équivalent au break d'avant). checkpoints = sorted((radii_steps[r_idx], r_idx) for r_idx in range(n_radii)) last_step = checkpoints[-1][0] # Process one direction at a time to limit GPU memory. # Store results as flat CPU arrays — transfer back to GPU at the end. pos_results = [None] * n_dirs neg_results = [None] * n_dirs for d_idx in range(n_dirs): ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] # Pre-compute valid steps for this direction (jusqu'au dernier checkpoint) valid_steps = [] for step in range(1, last_step + 1): px = int(round(ddx * step)) py = int(round(ddy * step)) dist_m = math.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) if dist_m < res * 0.5: continue valid_steps.append((step, px, py, dist_m)) # Max cumulé des tangentes (float32 : moitié de VRAM vs float64) running_pos = xp.zeros((rows, cols), dtype=np.float32) running_neg = xp.zeros((rows, cols), dtype=np.float32) snapshots = {} cp_queue = list(checkpoints) for step, px, py, dist_m in valid_steps: # Slice from padded array, subtract original dem view = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] elev_diff = view - dem del view # free slice reference # Tangentes des angles (positive : terrain au-dessus, négative : en # dessous). fmax propage le non-NaN : le bord de padding ne compte # pas, comme avec l'ancien where(isnan) — en une seule opération. running_pos = xp.fmax(running_pos, xp.maximum(elev_diff, 0) / dist_m) running_neg = xp.fmax(running_neg, xp.maximum(-elev_diff, 0) / dist_m) del elev_diff # free intermediate # Snapshot du checkpoint atteint : conversion en angle UNE fois while cp_queue and step >= cp_queue[0][0]: _, r_idx = cp_queue.pop(0) snapshots[r_idx] = (xp.arctan(running_pos), xp.arctan(running_neg)) if not cp_queue: break # Checkpoints jamais atteints (steps invalides) : état final du balayage while cp_queue: _, r_idx = cp_queue.pop(0) snapshots[r_idx] = (xp.arctan(running_pos), xp.arctan(running_neg)) # Store results on CPU, free GPU memory before next direction pos_results[d_idx] = to_cpu(xp.stack([snapshots[r][0] for r in range(n_radii)])) neg_results[d_idx] = to_cpu(xp.stack([snapshots[r][1] for r in range(n_radii)])) del running_pos, running_neg, snapshots gpu_cleanup() # Free the large padded array del padded gpu_cleanup() # Reassemble into final arrays (on CPU to avoid GPU memory pressure) pos_angles = np.array(pos_results) neg_angles = np.array(neg_results) return pos_angles, neg_angles def _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m=None): """Ray-tracing avec repli CPU automatique si la VRAM est insuffisante. Les dalles 0,2 m (5000×5000 px) multi-rayons peuvent dépasser la VRAM disponible (GPU partagé avec d'autres services) : plutôt que d'abandonner la visualisation, on désactive le GPU pour ce worker et on relance le calcul sur CPU. """ try: return _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m) except Exception as e: if _gpu_mod.is_gpu_active() and "out of memory" in str(e).lower(): logger.warning(" ⚠ VRAM insuffisante (ray-tracing) — repli CPU pour ce worker") _gpu_mod.disable_gpu() gpu_cleanup() return _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m) raise # ============================================================ # Core terrain visualizations # ============================================================ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None): """Generate multi-directional hillshade with contrast enhancement — GPU if available. Combines 8-direction hillshade with slope shading for balanced illumination. Applies percentile normalization and gamma correction to restore contrast lost by averaging multiple azimuths. """ 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" used_gpu = _gpu_mod.HAS_GPU try: if shared: transform = shared.transform crs = shared.crs # Pas de copie GPU du DEM brut ici : seuls gradient/pente/aspect # servent (déjà partagés) — ~100 Mo de VRAM économisés par dalle 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) else: dem_np, transform, crs = _read_dem(dem_file) dem = to_gpu(dem_np) # Espacement = résolution (m/px) : sans lui la pente est fausse dy, dx = xp.gradient(dem, float(resolution) if resolution else 1.0) slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) aspect = xp.arctan2(dy, dx) sin_slope = xp.sin(slope) cos_slope = xp.cos(slope) # 8 azimuths for balanced illumination (eliminates directional bias) azimuts = [0, 45, 90, 135, 180, 225, 270, 315] altitude = 35 # Higher altitude for better micro-relief detection hillshades = [] alt_rad = xp.radians(xp.array(altitude)) sin_alt = xp.sin(alt_rad) cos_alt = xp.cos(alt_rad) for az in azimuts: az_rad = xp.radians(xp.array(az)) hs = sin_alt * sin_slope + cos_alt * cos_slope * xp.cos(az_rad - aspect) hillshades.append(xp.clip(hs, 0, 1)) combined_hillshade = xp.mean(xp.array(hillshades), axis=0) slope_shaded = cos_slope combined = 0.7 * combined_hillshade + 0.3 * slope_shaded # Contrast enhancement: percentile stretch + gamma combined_np = to_cpu(combined) 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) if p98 - p2 > 0.01: combined_np = np.clip((combined_np - p2) / (p98 - p2), 0, 1) # Gamma correction to enhance shadows gamma = 0.8 combined_np = np.power(combined_np, gamma) _save_tif(output, combined_np.astype(np.float32), transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur hillshade: {e}", exc_info=True) return None def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): """Generate slope map (degrees) — GPU if available.""" 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" try: if shared: transform = shared.transform crs = shared.crs slope = shared.slope_deg nan_mask = shared.nan_mask if _gpu_mod.HAS_GPU: slope = to_gpu(slope) else: dem_np, transform, crs = _read_dem(dem_file) dem = to_gpu(dem_np) # Espacement = résolution (m/px) : sans lui la pente est fausse dy, dx = xp.gradient(dem, float(resolution) if resolution else 1.0) 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 _gpu_mod.HAS_GPU else slope, transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur slope: {e}", exc_info=True) return None def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): """Generate aspect (slope orientation) map — GPU if available. 0° = North, 90° = East, 180° = South, 270° = West. Direction toward which the terrain descends. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Aspect (Orientation des pentes){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_aspect.tif" try: if shared: transform = shared.transform crs = shared.crs dy = shared.dy dx = shared.dx nan_mask = shared.nan_mask if _gpu_mod.HAS_GPU: dy = to_gpu(dy) dx = to_gpu(dx) else: dem_np, transform, crs = _read_dem(dem_file) dem = to_gpu(dem_np) # Espacement = résolution (m/px) : sans lui la pente est fausse dy, dx = xp.gradient(dem, float(resolution) if resolution else 1.0) nan_mask = np.isnan(dem_np) aspect = xp.arctan2(dy, dx) * 180 / xp.pi aspect = xp.mod(aspect, 360) _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]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur aspect: {e}", exc_info=True) return None # ============================================================ # GPU-accelerated visualizations # ============================================================ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): """Sky-View Factor - ray-tracing on 16 azimuths, multi-radius (GPU if available). Traces rays in 16 directions at 3 radii (25, 50, 100m) and combines with weights favoring medium range for archaeological feature detection. SVF = (1/N) * sum(cos²(horizon_angle)). Valleys/crevices have low SVF, ridges/peaks have high SVF. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Sky-View Factor (ray-tracing multi-rayon){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_svf.tif" try: dem, dem_np, rows, cols, res, nan_mask, transform, crs = \ _prepare_dem_for_raycast(dem_file, shared, resolution) radii_m = [25, 50, 100] radius_weights = [0.3, 0.4, 0.3] # Medium range weighted more max_dist = int(max(radii_m) / res) # rayon réel en pixels, non tronqué n_dirs = 16 pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m) # pos/neg are now numpy arrays (CPU) — combine on CPU svf_combined = np.zeros((rows, cols), dtype=np.float32) for r_idx in range(len(radii_m)): horizon = np.maximum(pos_angles[:, r_idx], neg_angles[:, r_idx]) svf_r = np.mean(np.cos(horizon) ** 2, axis=0) svf_combined += svf_r * radius_weights[r_idx] svf_np = svf_combined svf_np[nan_mask] = np.nan _save_tif(output, svf_np, transform, crs) logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur SVF: {e}", exc_info=True) return None # Sous-échantillonnage du calcul d'openness : le lancé de rayons est l'étape # la plus coûteuse du pipeline (rayons de 100 m = 500 pas à 0,2 m/px). L'openness # étant un champ angulaire lisse (moyenne d'horizons jusqu'à 100 m), on la # calcule sur une grille décimée par blocs — max local pour l'openness # positive, min pour la négative, afin de préserver les reliefs qui bornent # l'angle d'horizon (talus, murs) — puis on la rééchantillonne à la # résolution demandée. Coût ÷ facteur³ (cellules ÷ facteur², rayons ÷ facteur) # ; facteur 2 ≈ ×8 plus rapide, rendu quasi identique. 1 = pleine résolution. OPENNESS_DOWNSAMPLE = 2 # Références de normalisation de l'openness (degrés : moyenne, écart-type). # Le z-score par tuile rendait l'échelle non jointive — même ouverture # physique, couleur différente d'une dalle à l'autre selon le relief # environnant. Médianes inter-tuiles mesurées à 0,2 m (décimation ×2) sur # des dalles réparties sur le territoire (plaine, bocage, forêt, montagne, # volcans, delta, urbain, littoral). Références FIGÉES : même ouverture = # même couleur sur toutes les tuiles. # Rayons du lancé de rayons (mètres), moyennés à poids égaux. OPENNESS_RADII_M = (25, 50, 100) OPENNESS_POS_REF = (5.571, 4.112) OPENNESS_NEG_REF = (5.974, 4.070) def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None): """Positive/Negative Openness - multi-radius ray-tracing with std normalization. Traces rays in 8 directions at 3 radii (25, 50, 100m) on a block-decimated grid (cf. OPENNESS_DOWNSAMPLE), then bilinearly resamples the result back to the requested resolution. Results are combined with equal weight across radii, then normalized by FIXED references (OPENNESS_POS_REF / OPENNESS_NEG_REF, degrees) so that adjacent tiles share one colour scale. """ name = "positive_openness" if positive else "negative_openness" gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing multi-rayon){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_{name}.tif" try: dem, dem_np, rows, cols, res, nan_mask, transform, crs = \ _prepare_dem_for_raycast(dem_file, shared, resolution) radii_m = list(OPENNESS_RADII_M) n_dirs = 8 full_rows, full_cols = rows, cols # Grille décimée par blocs (cf. OPENNESS_DOWNSAMPLE) factor = max(1, int(OPENNESS_DOWNSAMPLE)) if factor > 1 and rows >= 2 * factor and cols >= 2 * factor: r2 = (rows // factor) * factor c2 = (cols // factor) * factor blocks = dem[:r2, :c2].reshape(r2 // factor, factor, c2 // factor, factor) dem = blocks.max(axis=(1, 3)) if positive else blocks.min(axis=(1, 3)) rows, cols = dem.shape res = res * factor logger.info(f" Grille décimée ×{factor} ({rows}×{cols}) — rééchantillonnage final") max_dist = int(max(radii_m) / res) # rayon réel en pixels, non tronqué pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m) # Select positive or negative if positive: angles = pos_angles else: angles = neg_angles # Mean across directions and radii (equal weight) — on CPU now openness = np.mean(angles, axis=(0, 1)) openness_result = np.degrees(openness).astype(np.float32) # Retour à la grille pleine résolution (bilinéaire ; bord ajusté au plus # proche voisin si les dimensions n'étaient pas divisibles par le facteur) if (rows, cols) != (full_rows, full_cols): from scipy.ndimage import zoom zoomed = zoom(openness_result, factor, order=1) if zoomed.shape != (full_rows, full_cols): zoomed = np.pad(zoomed, ((0, full_rows - zoomed.shape[0]), (0, full_cols - zoomed.shape[1])), mode='edge') openness_result = zoomed.astype(np.float32) openness_result[nan_mask] = np.nan # Écart aux références figées, en sigmas : même ouverture physique = # même valeur sur toutes les tuiles → mosaïque de couleur jointive ref_mean, ref_std = OPENNESS_POS_REF if positive else OPENNESS_NEG_REF openness_result = (openness_result - ref_mean) / ref_std _save_tif(output, openness_result, transform, crs) logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur openness: {e}", exc_info=True) return None # ============================================================ # Relief orienté : openness locale × aspect en une seule image RGB # ============================================================ # # Clarté (CIELAB L*) = micro-relief : openness positive calculée sur le MNT # détendancé (MNT − gaussienne 10 m), rayons courts 5/10/20 m, plus un léger # ombrage directionnel. Teinte = orientation de la pente (aspect) sur le # cercle CIELAB, à clarté constante : aucune couleur ne crée de faux relief. # Pas de statistique par dalle (échelle log FIXE) et un support total # (rayon 20 m + lissage 4σ = 40 m) inférieur à la bande de raccord de 100 m : # les dalles adjacentes se raccordent sans couture. # # Coût maîtrisé : détendance et lancé de rayons tournent sur une grille # décimée à ~0,8 m (bloc max, cf. OPENNESS_DOWNSAMPLE) — seule l'openness # finale est agrandie ; noyau dédié qui n'accumule que la moyenne des angles # (CuPy RawKernel sur GPU, numba parallèle sur CPU, numpy en dernier # recours) ; couleurs par table pré-calculée (L* × teinte) au lieu d'une # conversion Lab → sRGB par pixel. RELIEF_DETREND_M = 10.0 # σ de la gaussienne retirée au MNT RELIEF_RADII_M = (5, 10, 20) # rayons du lancé de rayons, poids égaux RELIEF_N_DIRS = 16 # 16 directions : pas de losanges à 8 branches RELIEF_GRID_M = 0.8 # pas de la grille de calcul décimée RELIEF_OPEN_RANGE = (0.5, 12.0) # degrés, échelle log fixe (jointive) RELIEF_CHROMA = 60.0 # chroma CIELAB max (atteint vers L* = 50) RELIEF_SHADE_WEIGHT = 0.35 # part de l'ombrage directionnel dans L* # Azimut de l'ombrage, exprimé dans le repère de l'aspect ci-dessous # (arctan2(dy, dx), dy vers le sud) — réglage validé sur la galerie de rendus. RELIEF_SHADE_AZIMUTH = 315.0 RELIEF_SHADE_ALTITUDE = 45.0 RELIEF_NODATA_RGB = (38, 38, 41) _RELIEF_LUT = None def _lab_to_srgb(L, a, b): """CIELAB (D65) → sRGB [0, 1], écrêté.""" fy = (L + 16) / 116 fx = fy + a / 500 fz = fy - b / 200 finv = lambda t: np.where(t > 6 / 29, t ** 3, 3 * (6 / 29) ** 2 * (t - 4 / 29)) X, Y, Z = 0.95047 * finv(fx), finv(fy), 1.08883 * finv(fz) rgb = np.stack([3.2406 * X - 1.5372 * Y - 0.4986 * Z, -0.9689 * X + 1.8758 * Y + 0.0415 * Z, 0.0557 * X - 0.2040 * Y + 1.0570 * Z], -1) rgb = np.where(rgb <= 0.0031308, 12.92 * rgb, 1.055 * np.clip(rgb, 0, None) ** (1 / 2.4) - 0.055) return np.clip(rgb, 0, 1) def _relief_lut(): """Table (256 niveaux de L*, 360 teintes) → sRGB uint8, calculée une fois. La chroma décroît vers le noir et le blanc (L*(100−L*)/2500) : les couleurs restent dans la gamme sRGB au lieu d'être écrêtées. """ global _RELIEF_LUT if _RELIEF_LUT is None: L = np.linspace(0, 100, 256)[:, None] h = np.radians(np.arange(360))[None, :] C = RELIEF_CHROMA * np.clip(L * (100 - L) / 2500.0, 0, 1) _RELIEF_LUT = (_lab_to_srgb(np.broadcast_to(L, (256, 360)), C * np.cos(h), C * np.sin(h)) * 255 + 0.5).astype(np.uint8) return _RELIEF_LUT def _horizon_rays(res, n_dirs, radii_m): """Décalages (col, ligne) de chaque pas de rayon, distances et checkpoints. Même géométrie que _ray_trace_horizons_core : direction k à l'angle 2πk/n, pas arrondis au pixel, distance = pas × résolution. """ max_step = max(1, int(max(radii_m) / res)) angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) steps = np.arange(1, max_step + 1) offs = np.empty((n_dirs, max_step, 2), dtype=np.int32) offs[:, :, 0] = np.rint(np.cos(angles)[:, None] * steps) offs[:, :, 1] = np.rint(np.sin(angles)[:, None] * steps) dist = (steps * res).astype(np.float32) cps = np.array(sorted(max(1, min(int(r / res), max_step)) for r in radii_m), dtype=np.int32) return offs, dist, cps _HORIZON_CUDA_SRC = r''' extern "C" __global__ void mean_horizon(const float* dem, const int* offs, const float* dist, const int* cps, const int rows, const int cols, const int nd, const int ns, const int nr, float* out) { long idx = (long)blockDim.x * blockIdx.x + threadIdx.x; if (idx >= (long)rows * cols) return; int i = idx / cols, j = idx % cols; float z0 = dem[idx], acc = 0.f; for (int d = 0; d < nd; ++d) { float run = 0.f; int k = 0; for (int s = 0; s < ns; ++s) { int ii = i + offs[(d * ns + s) * 2 + 1]; int jj = j + offs[(d * ns + s) * 2]; if (ii >= 0 && ii < rows && jj >= 0 && jj < cols) run = fmaxf(run, (dem[(long)ii * cols + jj] - z0) / dist[s]); while (k < nr && s + 1 >= cps[k]) { acc += atanf(run); ++k; } } while (k < nr) { acc += atanf(run); ++k; } } out[idx] = acc / (float)(nd * nr); } ''' _horizon_cuda_kernel = None _horizon_numba_kernel = None def _mean_horizon_gpu(dem, offs, dist, cps): """Un thread CUDA par pixel ; tout reste en registres (aucun tableau intermédiaire par direction ou par rayon).""" global _horizon_cuda_kernel cp = _gpu_mod._cp if _horizon_cuda_kernel is None: _horizon_cuda_kernel = cp.RawKernel(_HORIZON_CUDA_SRC, 'mean_horizon') rows, cols = dem.shape d_dem = cp.ascontiguousarray(cp.asarray(dem, dtype=cp.float32)) out = cp.empty((rows, cols), dtype=cp.float32) n = rows * cols threads = 256 _horizon_cuda_kernel(((n + threads - 1) // threads,), (threads,), (d_dem, cp.asarray(offs), cp.asarray(dist), cp.asarray(cps), np.int32(rows), np.int32(cols), np.int32(offs.shape[0]), np.int32(offs.shape[1]), np.int32(len(cps)), out)) return out def _mean_horizon_numba(dem, offs, dist, cps): """Même noyau en numba parallèle (une ligne de pixels par thread CPU).""" global _horizon_numba_kernel if _horizon_numba_kernel is None: from numba import njit, prange @njit(parallel=True, cache=True, fastmath=True) def _kernel(dem, offs, dist, cps): rows, cols = dem.shape nd, ns, nr = offs.shape[0], offs.shape[1], cps.shape[0] out = np.empty((rows, cols), dtype=np.float32) for i in prange(rows): for j in range(cols): z0 = dem[i, j] acc = 0.0 for d in range(nd): run = 0.0 k = 0 for s in range(ns): ii = i + offs[d, s, 1] jj = j + offs[d, s, 0] if ii >= 0 and ii < rows and jj >= 0 and jj < cols: t = (dem[ii, jj] - z0) / dist[s] if t > run: run = t while k < nr and s + 1 >= cps[k]: acc += math.atan(run) k += 1 while k < nr: acc += math.atan(run) k += 1 out[i, j] = acc / (nd * nr) return out _horizon_numba_kernel = _kernel return _horizon_numba_kernel(np.ascontiguousarray(dem, dtype=np.float32), offs, dist, cps) def _mean_horizon_numpy(dem, offs, dist, cps): """Repli vectorisé (sans numba ni GPU) : un décalage de tableau par pas.""" rows, cols = dem.shape pad = int(np.abs(offs).max()) padded = np.pad(dem.astype(np.float32), pad, constant_values=np.nan) acc = np.zeros((rows, cols), dtype=np.float64) for d in range(offs.shape[0]): run = np.zeros((rows, cols), dtype=np.float32) k = 0 for s in range(offs.shape[1]): px, py = offs[d, s] view = padded[pad + py:pad + py + rows, pad + px:pad + px + cols] run = np.fmax(run, (view - dem) / dist[s]) while k < len(cps) and s + 1 >= cps[k]: acc += np.arctan(run) k += 1 while k < len(cps): acc += np.arctan(run) k += 1 return (acc / (offs.shape[0] * len(cps))).astype(np.float32) def _mean_horizon_angle(dem, res, n_dirs, radii_m): """Moyenne (directions × rayons) de l'angle d'horizon positif, en radians. GPU (CuPy RawKernel) si disponible — le résultat reste alors sur le GPU —, sinon numba, sinon numpy. Un échec GPU (compilation NVRTC, VRAM) bascule sur le CPU pour ce calcul sans perdre la visualisation. """ offs, dist, cps = _horizon_rays(res, n_dirs, radii_m) if _gpu_mod.is_gpu_active(): try: return _mean_horizon_gpu(dem, offs, dist, cps), "GPU" except Exception as e: logger.warning(f" ⚠ Noyau GPU relief orienté indisponible ({e}) — repli CPU") dem = to_cpu(dem) try: return _mean_horizon_numba(dem, offs, dist, cps), "numba" except ImportError: return _mean_horizon_numpy(dem, offs, dist, cps), "numpy" def _pad_to(m, arr, rows, cols): """Complète par répétition du bord (dimensions non multiples du facteur).""" if arr.shape == (rows, cols): return arr return m.pad(arr[:rows, :cols], ((0, max(0, rows - arr.shape[0])), (0, max(0, cols - arr.shape[1]))), mode='edge') def _relief_params(): """Constantes scalaires du rendu (ombrage sans trigonométrie par pixel). cos(pente) = 1/√(1+g²) et sin(pente)·cos(az − aspect) = (cos az·dx + sin az·dy)/√(1+g²), avec aspect = arctan2(dy, dx) : l'ombrage ne coûte qu'une racine par pixel. """ zen = math.radians(90.0 - RELIEF_SHADE_ALTITUDE) az = math.radians(RELIEF_SHADE_AZIMUTH) lo, hi = RELIEF_OPEN_RANGE return (math.cos(zen), math.sin(zen) * math.cos(az), math.sin(zen) * math.sin(az), lo, 1.0 / math.log(hi / lo), RELIEF_SHADE_WEIGHT) def _relief_colorize_xp(m, openness, dx, dy, lut): """Colorisation vectorisée (CuPy sur GPU, numpy en repli).""" cz, sa_x, sa_y, lo, inv_log, w = _relief_params() inv = 1.0 / m.sqrt(1.0 + dx * dx + dy * dy) shade = m.clip((cz + sa_x * dx + sa_y * dy) * inv / cz, 0, 1.25) / 1.25 t = m.clip(m.log(m.maximum(openness, 1e-3) / lo) * inv_log, 0, 1) L = 12.0 + 84.0 * ((1 - w) * (1 - t) + w * shade) del inv, shade, t li = m.clip(m.rint(L * 2.55), 0, 255).astype(m.int32) hue = m.mod(m.rint(m.degrees(m.arctan2(dy, dx))), 360).astype(m.int32) return lut[li, hue] _relief_color_numba_kernel = None def _relief_colorize_numba(open_c, f, dx, dy, lut): """Colorisation CPU en une passe (numba parallèle) : l'openness est lue sur la grille décimée par interpolation bilinéaire alignée sur les centres de pixels (équivalent de xp_zoom), sans tableau pleine résolution intermédiaire.""" global _relief_color_numba_kernel if _relief_color_numba_kernel is None: from numba import njit, prange @njit(parallel=True, cache=True, fastmath=True) def _kernel(open_c, f, dx, dy, lut, cz, sa_x, sa_y, lo, inv_log, w): rows, cols = dx.shape rc, cc = open_c.shape out = np.empty((rows, cols, 3), dtype=np.uint8) for i in prange(rows): u = (i + 0.5) / f - 0.5 u = min(max(u, 0.0), rc - 1.0) i0 = int(u) i1 = min(i0 + 1, rc - 1) fu = u - i0 for j in range(cols): v = (j + 0.5) / f - 0.5 v = min(max(v, 0.0), cc - 1.0) j0 = int(v) j1 = min(j0 + 1, cc - 1) fv = v - j0 o = ((1 - fu) * ((1 - fv) * open_c[i0, j0] + fv * open_c[i0, j1]) + fu * ((1 - fv) * open_c[i1, j0] + fv * open_c[i1, j1])) gx = dx[i, j] gy = dy[i, j] inv = 1.0 / math.sqrt(1.0 + gx * gx + gy * gy) sh = (cz + sa_x * gx + sa_y * gy) * inv / cz sh = min(max(sh, 0.0), 1.25) / 1.25 t = math.log(max(o, 1e-3) / lo) * inv_log t = min(max(t, 0.0), 1.0) L = 12.0 + 84.0 * ((1 - w) * (1 - t) + w * sh) li = min(max(int(round(L * 2.55)), 0), 255) h = int(round(math.degrees(math.atan2(gy, gx)))) % 360 out[i, j, 0] = lut[li, h, 0] out[i, j, 1] = lut[li, h, 1] out[i, j, 2] = lut[li, h, 2] return out _relief_color_numba_kernel = _kernel return _relief_color_numba_kernel(np.ascontiguousarray(open_c), float(f), np.ascontiguousarray(dx, dtype=np.float32), np.ascontiguousarray(dy, dtype=np.float32), lut, *_relief_params()) def generate_relief_oriente(dem_file, basename, vis_dir, resolution, shared=None): """Relief orienté : image RGB unique fusionnant openness locale et aspect. L* (clarté) = openness positive locale (MNT détendancé, rayons courts) 65 % + ombrage directionnel 35 % ; teinte = orientation de la pente ; chroma CIELAB fixe. Sortie : GeoTIFF RGB uint8 (rendu tel quel, comme les fonds IGN). """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Relief orienté (openness locale × aspect){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_relief_oriente.tif" try: dem, dem_np, rows, cols, res, nan_mask, transform, crs = \ _prepare_dem_for_raycast(dem_file, shared, resolution) res = float(res) use_gpu = _gpu_mod.is_gpu_active() filled = (shared.filled_gpu if shared is not None and use_gpu else None) if filled is None: filled = to_gpu(dem) if use_gpu else dem.astype(np.float32) # 1. Grille décimée : bloc max (préserve les reliefs qui bornent # l'horizon) et bloc moyen (support de la tendance à retirer) f = max(1, int(round(RELIEF_GRID_M / res))) if rows < 2 * f or cols < 2 * f: f = 1 r2, c2 = (rows // f) * f, (cols // f) * f blocks = filled[:r2, :c2].reshape(r2 // f, f, c2 // f, f) coarse_max = blocks.max(axis=(1, 3)) coarse_mean = blocks.mean(axis=(1, 3)) res_c = res * f # 2. Détendance (MNT − gaussienne) sur la grille décimée trend = xp_gaussian_filter(coarse_mean, RELIEF_DETREND_M / res_c) detrended = coarse_max - trend del blocks, coarse_max, coarse_mean, trend t_prep = time.time() - t0 # 3. Openness locale : angle d'horizon moyen (radians) sur la grille décimée t1 = time.time() angle, engine = _mean_horizon_angle(detrended, res_c, RELIEF_N_DIRS, RELIEF_RADII_M) t_ray = time.time() - t1 on_gpu = engine == "GPU" del detrended # 4. Couleurs : agrandissement de l'openness, ombrage, table L* × teinte t2 = time.time() if shared is not None: dy, dx = shared.dy, shared.dx else: dy, dx = np.gradient(to_cpu(filled).astype(np.float32), res) lut = _relief_lut() if on_gpu: openness = _gpu_mod.xp_zoom(xp.degrees(angle), f, order=1) if f > 1 else xp.degrees(angle) openness = _pad_to(xp, xp.asarray(openness), rows, cols) rgb = to_cpu(_relief_colorize_xp(xp, openness, to_gpu(dx), to_gpu(dy), xp.asarray(lut))) del openness gpu_cleanup() engine_color = "GPU" else: angle = to_cpu(angle) try: rgb = _relief_colorize_numba(np.degrees(angle).astype(np.float32), f, dx, dy, lut) engine_color = "numba" except ImportError: openness = _gpu_mod.xp_zoom(np.degrees(angle), f, order=1) if f > 1 else np.degrees(angle) rgb = _relief_colorize_xp(np, _pad_to(np, openness, rows, cols), dx, dy, lut) engine_color = "numpy" del angle, dx, dy t_color = time.time() - t2 rgb[nan_mask] = RELIEF_NODATA_RGB _save_tif(output, np.moveaxis(rgb, -1, 0), transform, crs, dtype='uint8', count=3) logger.info(f" ✓ Relief orienté terminé ({time.time()-t0:.1f}s : " f"préparation {t_prep:.1f}s, rayons {t_ray:.1f}s [{engine}], " f"couleurs {t_color:.1f}s [{engine_color}])") return output except Exception as e: logger.error(f" ✗ Erreur relief orienté: {e}", exc_info=True) return None # Densité de points sol : 16 niveaux sur une échelle log fixe (même gris = # même densité partout, mosaïque jointive). Le niveau k commence à # DENSITY_BASE × 2^(k/2) pts/m² : 0 = ≤ 0,35 (noir, y compris sans aucun # point), 15 = ≥ 45 (blanc) ; deux niveaux = densité doublée. DENSITY_LEVELS = 16 DENSITY_BASE = 0.25 def density_levels(density): """Niveau 0..15 d'une densité (pts/m²).""" with np.errstate(divide="ignore", invalid="ignore"): k = np.floor(2.0 * np.log2(np.asarray(density, dtype=np.float64) / DENSITY_BASE)) return np.clip(np.nan_to_num(k, nan=0.0, neginf=0.0), 0, DENSITY_LEVELS - 1).astype(np.uint8) def generate_densite_sol(dem_file, basename, vis_dir, resolution, shared=None): """Densité des points sol retenus pour le MNT, en 16 niveaux de gris. Lue dans le fichier annexe écrit à la rastérisation (dtm.density_path) : grille de 1 m, gardée telle quelle (la densité n'a pas plus de détail — image 25× plus légère qu'à 0,2 m). Sortie : niveau 0..15 (float32). """ logger.info(" → Densité de points sol...") t0 = time.time() output = vis_dir / f"{basename}_densite_sol.tif" try: from .dtm import density_path src_path = density_path(dem_file) if not src_path.exists(): logger.error(f" ✗ Densité absente ({src_path.name}) : régénérer le DTM") return None with rasterio.open(src_path) as src: density = src.read(1) transform, crs = src.transform, src.crs _save_tif(output, density_levels(density).astype(np.float32), transform, crs) logger.info(f" ✓ Densité de points terminée ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Erreur densité de points: {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 _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" 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) # Adaptive scales: finer at higher resolution min_scale = max(2.0, resolution * 4) # Archaeological scales: focus on 2-50m range. # Small features (ditches, walls, post-holes) need 2-10m. # Medium features (enclosures, roundhouses) need 10-25m. # Large scales (50m+) are kept only for context with low weight. candidate_scales = [2, 3, 5, 8, 10, 15, 25, 50] sigmas = [s for s in candidate_scales if s >= min_scale] # Weights: favor small-to-medium scales where archaeo features live scale_weights = { 2: 1.5, 3: 1.8, 5: 2.0, 8: 1.8, 10: 1.5, 15: 1.3, 25: 1.0, 50: 0.5, } weights = np.array([scale_weights.get(s, 1.0) for s in sigmas]) logger.info(f" MSRM échelles: {sigmas}m") lrm_stack = [] for sigma in sigmas: 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 valid_lrm = lrm[~nan_mask] lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 lrm = lrm / lrm_std # Clip |z| to 3.0 to prevent large-scale outliers from drowning small features lrm = np.clip(lrm, -3.0, 3.0) lrm_stack.append(lrm.astype(np.float32)) # Weighted combination — preserve sign for RdBu_r colormap # Positive = elevated (red), Negative = depression (blue) 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') # Signed RMS: magnitude from RMS, sign from weighted mean signed_mean = np.nansum(lrm_array * weights_3d, axis=0) / np.sum(weights) rms_magnitude = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights)) mslrm = np.sign(signed_mean) * rms_magnitude mslrm[nan_mask] = np.nan _save_tif(output, mslrm.astype(np.float32), transform, crs) logger.info(f" ✓ MSRM terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur MSRM: {e}", exc_info=True) return None # ============================================================ # SAILORE # ============================================================ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None): """SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available). Kernel size adapts to local slope: flat areas get larger kernels, steep areas get smaller kernels. Scales adapt to resolution. Reuses shared.lrm_15 when available to avoid recomputation. """ 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" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np nan_mask = shared.nan_mask slope_deg = shared.slope_deg else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) gy, gx = np.gradient(dem_np, resolution) slope = np.arctan(np.sqrt(gx**2 + gy**2)) slope_deg = np.degrees(slope) slope_deg[nan_mask] = np.nan # 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 slope_norm = np.clip(slope_deg / 30.0, 0, 1) # LRM fine (2m) — always compute if shared: lrm_fine = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_min) else: lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min) lrm_fine[nan_mask] = np.nan # LRM medium (13.5m) — reuse shared.lrm_15 (σ=15m) when available sigma_mid = (sigma_min + sigma_max) / 2 if shared and abs(15.0 / resolution - sigma_mid) < 2.0 / resolution: # shared.lrm_15 is close enough to medium scale lrm_medium = shared.lrm_15.copy() else: if shared: lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_mid) else: lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_mid) lrm_medium[nan_mask] = np.nan # LRM coarse (25m) — always compute if shared: lrm_coarse = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_max) else: lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max) lrm_coarse[nan_mask] = np.nan w_fine = slope_norm w_medium = 1 - 2 * np.abs(slope_norm - 0.5) w_coarse = 1 - slope_norm w_total = w_fine + w_medium + w_coarse w_total[w_total == 0] = 1 sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total sailore[nan_mask] = np.nan # Z-score (σ locales) : unités comparables entre tuiles, plage de # rendu fixe ±3σ → mosaïque de couleur homogène valid = sailore[~nan_mask] if len(valid) > 0: std_val = max(np.nanstd(valid), 0.01) sailore = (sailore - np.nanmean(valid)) / std_val _save_tif(output, sailore.astype(np.float32), transform, crs) logger.info(f" ✓ SAILORE terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur SAILORE: {e}", exc_info=True) return None # ============================================================ # Roughness # ============================================================ def _integral_sums(x): """Sommes intégrales 2D : S[i,j] = somme de x[0:i, 0:j] (float64). Ligne/colonne 0 remplies de zéros — permet la somme d'une fenêtre quelconque par 4 coins, y compris contre le bord (indices 0). """ rows, cols = x.shape S = xp.zeros((rows + 1, cols + 1), dtype=np.float64) S[1:, 1:] = xp.cumsum(xp.cumsum(x.astype(np.float64), axis=0), axis=1) return S def _box_std_from_integral(Sx, Sx2, size): """Écart-type local sur fenêtre size×size via sommes intégrales. Coût indépendant de la taille de fenêtre (4 accès par pixel) — contre un uniform_filter dont le coût croît avec la fenêtre (75 px à 0,2 m pour l'échelle large). Aux bords, la fenêtre est tronquée et normalisée par le nombre réel d'éléments (les dalles se recouvrent, le bord est sans effet visuel). """ rows = Sx.shape[0] - 1 cols = Sx.shape[1] - 1 r = size // 2 iy0 = xp.maximum(xp.arange(rows) - r, 0) iy1 = xp.minimum(xp.arange(rows) + r + 1, rows) ix0 = xp.maximum(xp.arange(cols) - r, 0) ix1 = xp.minimum(xp.arange(cols) + r + 1, cols) # Nombre d'éléments réels de la fenêtre (tronquée aux bords) counts = ((iy1 - iy0)[:, None] * (ix1 - ix0)[None, :]).astype(np.float64) def box_sum(S): return (S[iy1][:, ix1] - S[iy0][:, ix1] - S[iy1][:, ix0] + S[iy0][:, ix0]) mean = box_sum(Sx) / counts mean_sq = box_sum(Sx2) / counts return xp.sqrt(xp.maximum(mean_sq - mean * mean, 0)) # Références de normalisation de la rugosité (mètres d'écart-type local). # Médianes inter-tuiles mesurées sur 20 dalles réelles à 0,2 m : la # normalisation par tuile (z-score) rendait l'échelle non jointive — l'écart # variait de 0,05 à 0,58 m selon la tuile pour l'échelle fine. References # FIGÉES : même rugosité physique = même valeur sur toutes les tuiles. ROUGHNESS_FINE_REF_M = 0.156 ROUGHNESS_BROAD_REF_M = 0.475 def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): """Surface roughness - multi-scale standard deviation (GPU-accelerated). Combines fine (3m) and broad (15m) roughness for better detection of archaeological features at multiple scales. Les écarts-types locaux sont calculés par sommes intégrales : deux cumsum partagés entre les deux échelles, extraction par 4 coins — coût constant quelle que soit la fenêtre (un uniform_filter coûte proportionnellement à sa taille, 75 px à 0,2 m pour l'échelle large). Normalisation par références FIXÉES (médianes mesurées) et non par tuile : les mosaïques sont jointives, même valeur = même couleur. """ 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" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np nan_mask = shared.nan_mask if _gpu_mod.HAS_GPU: filled = shared.filled_gpu else: filled = shared.filled else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) if _gpu_mod.HAS_GPU: filled = to_gpu(filled) # Sommes intégrales partagées par les deux échelles (X et X²) Sx = _integral_sums(filled) Sx2 = _integral_sums(filled.astype(np.float64) ** 2) fine_size = max(3, int(3 / resolution)) if fine_size % 2 == 0: fine_size += 1 broad_size = max(3, int(15 / resolution)) if broad_size % 2 == 0: broad_size += 1 roughness_fine = to_cpu(_box_std_from_integral(Sx, Sx2, fine_size)) roughness_broad = to_cpu(_box_std_from_integral(Sx, Sx2, broad_size)) del Sx, Sx2 gpu_cleanup() roughness_fine[nan_mask] = np.nan roughness_broad[nan_mask] = np.nan # Combinaison pondérée, échelle physique commune (références fixées : # jointive entre tuiles — cf. constantes module) roughness = (0.7 * roughness_fine / ROUGHNESS_FINE_REF_M + 0.3 * roughness_broad / ROUGHNESS_BROAD_REF_M) roughness[nan_mask] = np.nan _save_tif(output, roughness, transform, crs) logger.info(f" ✓ Rugosité terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur rugosité: {e}", exc_info=True) return None # ============================================================ # Exposition des surfaces (Éclairage Solaire) # ============================================================ def generate_solar(dem_file, basename, vis_dir, resolution, shared=None): """Generate solar irradiance simulation. Simulates morning sunlight (azimuth 90°, altitude 30°) to reveal subtle topographic features through shadow effects. """ logger.info(" → Exposition des surfaces (Éclairage Solaire)...") t0 = time.time() output = vis_dir / f"{basename}_solar.tif" try: if shared: transform = shared.transform crs = shared.crs nan_mask = shared.nan_mask dx = shared.dx dy = shared.dy else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) dem_filled, _ = _fill_nans(dem_np) dy, dx = np.gradient(dem_filled, resolution, resolution) # Solar parameters: morning sun (azimuth 90° = east, altitude 30°) sun_azimuth = np.radians(90) sun_altitude = np.radians(30) # Aspect from gradient aspect = np.degrees(np.arctan2(-dx, -dy)) aspect[aspect < 0] += 360 # Slope in radians slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) # Solar irradiance calculation irradiance = (np.sin(sun_altitude) * np.sin(slope_rad) + np.cos(sun_altitude) * np.cos(slope_rad) * np.cos(np.radians(aspect) - sun_azimuth)) # Clip to valid range [0, 1] irradiance = np.clip(irradiance, 0, 1) irradiance[nan_mask] = np.nan _save_tif(output, irradiance.astype(np.float32), transform, crs) logger.info(f" ✓ Exposition des surfaces terminée ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Erreur exposition: {e}", exc_info=True) return None # ============================================================ # Wavelet (Mexican Hat) # ============================================================ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): """Mexican Hat wavelet multi-scale analysis (GPU if available). Focused on small archaeological structures (paths, ditches, ramparts). CWT 2D at scales [1, 2, 5, 10, 20, 50]m (0.5m added below 0.25m/px). The 100m scale was dropped: it mostly responds to landforms (hills, valleys), not to structures. Large landforms are removed first by subtracting a Gaussian local-mean trend (~35m). The residual is analyzed relative to its ~35m neighborhood, so a ditch on a hilltop or slope does not stand out more than the same ditch on flat ground (topographic-position independence). A Gaussian smoothing preserves locally planar slopes, so slopes are removed too. Uses robust per-scale normalization (MAD) and median-centered weighted RMS combination with emphasis on small scales (1-10m). The output is an index relative to the tile's own median level (≈ 1): comparable from tile to tile, so a single fixed color stretch works. """ 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" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np nan_mask = shared.nan_mask filled = shared.filled.astype(np.float64) else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np.astype(np.float64)) min_scale = max(resolution * 2, 1.0) # 100m retirée : elle répond surtout aux grandes formes du terrain, # pas aux structures. Accent sur 1-10m (petites structures). candidate_scales = [0.5, 1, 2, 5, 10, 20, 50] scales = [s for s in candidate_scales if s >= min_scale] scale_weights = { 0.5: 0.7, 1.0: 1.2, 2.0: 1.8, 5.0: 2.2, 10.0: 2.0, 20.0: 1.5, 50.0: 0.8, } weights = np.array([scale_weights.get(s, 1.0) for s in scales]) logger.info(f" Échelles CWT: {scales}m (résolution {resolution}m/px)") from scipy.ndimage import gaussian_laplace, gaussian_filter # Retrait des grands volumes (collines, vallées) : on soustrait une # moyenne locale gaussienne avant la CWT. Un lissage gaussien préserve # les pentes planes, donc le résidu est analysé par rapport à son # voisinage ~35m : un fossé en sommet ou en flanc de colline ne # ressort pas plus que le même fossé à plat. # Fraction conservée pour une structure gaussienne de largeur σ_f : # σ_t²/(σ_f²+σ_t²) → 10m : 92%, 20m : 75%, colline 150m : 5%. # Mesuré sur MNT synthétique bruité : le contraste des petites # structures est insensible à σ_t ; seul le fond sommet/plat varie # (1.71 sans détendage → 1.21 à 35m). detrend_sigma_m = 35.0 detrend_sigma_px = detrend_sigma_m / resolution if _gpu_mod.HAS_GPU: try: from cupyx.scipy.ndimage import gaussian_filter as gpu_gaussian_filter trend = to_cpu(gpu_gaussian_filter(to_gpu(filled), sigma=detrend_sigma_px)) except Exception: trend = gaussian_filter(filled, sigma=detrend_sigma_px) else: trend = gaussian_filter(filled, sigma=detrend_sigma_px) residual = filled - trend del trend logger.info(f" Retrait des grands volumes (tendance gaussienne {detrend_sigma_m:.0f}m)") wavelet_stack = [] for scale_m in scales: sigma_px = scale_m / resolution if _gpu_mod.HAS_GPU: try: from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace response = -gpu_gaussian_laplace(to_gpu(residual), sigma=sigma_px) response = to_cpu(response) except Exception: response = -gaussian_laplace(residual, sigma=sigma_px) else: response = -gaussian_laplace(residual, sigma=sigma_px) response[nan_mask] = np.nan valid = response[~nan_mask] # σ robuste (MAD) : le std classique est gonflé par les queues # (structures marquées, bords de tuile) et varie fortement d'une # tuile à l'autre — cause première des dominantes de couleur # par tuile sur la carte. mad = np.nanmedian(np.abs(valid - np.nanmedian(valid))) if len(valid) > 0 else 0.0 response = response / max(1.4826 * mad, 0.01) wavelet_stack.append(response) stack = np.array(wavelet_stack) weights_3d = weights[:, np.newaxis, np.newaxis] with np.errstate(invalid='ignore', divide='ignore'): with warnings.catch_warnings(): warnings.filterwarnings('ignore', message='Mean of empty slice') combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(weights)) combined[nan_mask] = np.nan # Recentrage par la médiane de la tuile : la RMS devient un indice # relatif au niveau moyen de la tuile (médiane = 1). La distribution # est alors comparable d'une tuile à l'autre — condition pour un # étirement couleur global fixe et homogène entre tuiles. finite = combined[np.isfinite(combined)] if finite.size: combined = combined / max(float(np.median(finite)), 0.01) combined[nan_mask] = np.nan _save_tif(output, combined.astype(np.float32), transform, crs) logger.info(f" ✓ Ondelette terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Erreur ondelette: {e}", exc_info=True) return None # ============================================================ # Flow Accumulation helpers (module-level for numba caching) # ============================================================ 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. NaN cells are treated as closed (never pushed into the heap). """ result = _priority_flood_numba(dem, nodata_mask) if result is not None: return result return _priority_flood_python(dem, nodata_mask) def _priority_flood_python(dem, nodata_mask): """Pure-Python fallback for _priority_flood (used when numba is unavailable).""" import heapq rows, cols = dem.shape filled = dem.copy() closed = nodata_mask.copy() open_queue = [] 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 closed[nr, nc] = True heapq.heappush(open_queue, (filled[nr, nc], nr, nc)) return filled def _priority_flood_numba(dem, nodata_mask): """JIT-compiled priority-flood via binary min-heap (~200x faster than Python). Returns None if numba is unavailable (caller falls back to Python). """ try: from numba import njit except ImportError: return None @njit(cache=True) def _flood(dem, nodata): rows, cols = dem.shape filled = dem.copy() flat = filled.ravel() closed = nodata.copy() n = rows * cols heap = np.empty(n, dtype=np.int64) heap_size = 0 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) for r in range(rows): for c in (0, cols - 1): if not closed[r, c]: heap[heap_size] = r * cols + c heap_size += 1 closed[r, c] = True for c in range(1, cols - 1): for r in (0, rows - 1): if not closed[r, c]: heap[heap_size] = r * cols + c heap_size += 1 closed[r, c] = True while heap_size > 0: cell = heap[0] elev = flat[cell] heap_size -= 1 if heap_size > 0: heap[0] = heap[heap_size] i = 0 while True: l = 2 * i + 1 r = 2 * i + 2 smallest = i if l < heap_size and flat[heap[l]] < flat[heap[smallest]]: smallest = l if r < heap_size and flat[heap[r]] < flat[heap[smallest]]: smallest = r if smallest == i: break heap[i], heap[smallest] = heap[smallest], heap[i] i = smallest r = cell // cols c = cell % cols for d in range(8): nr = r + dy8[d] nc = c + dx8[d] if 0 <= nr < rows and 0 <= nc < cols and not closed[nr, nc]: ncell = nr * cols + nc if flat[ncell] < elev: flat[ncell] = elev closed[nr, nc] = True heap[heap_size] = ncell heap_size += 1 child = heap_size - 1 while child > 0: parent = (child - 1) // 2 if flat[heap[child]] < flat[heap[parent]]: heap[child], heap[parent] = heap[parent], heap[child] child = parent else: break return filled return _flood(dem, nodata_mask) def _d8_accumulate_numba(dem_filled, flow_dir, nodata_mask, rows, cols): """JIT-compiled D8 flow accumulation (top-down via elevation sort). Uses numba for ~100x speedup over pure Python loop. Falls back to pure Python if numba is unavailable. """ try: from numba import njit @njit(cache=True) def _accumulate(dem, fdir, nodata, 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) for r in range(rows): for c in range(cols): if nodata[r, c]: flow_acc[r, c] = 0.0 # Sort cells by elevation descending (highest first) n_cells = rows * cols flat_dem = dem.ravel() sort_idx = np.argsort(-flat_dem) # Accumulate top-down (highest cell first) for i in range(n_cells): cell = sort_idx[i] r = cell // cols c = cell % cols if nodata[r, c]: continue d = fdir[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[nr, nc]: flow_acc[nr, nc] += flow_acc[r, c] return flow_acc return _accumulate(dem_filled, flow_dir, nodata_mask, rows, cols) except ImportError: return None # ============================================================ # Flow Accumulation # ============================================================ def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=None): """Flow Accumulation — priority-flood sink filling + D8 accumulation. Detects channels, ditches, and drainage paths by computing how many upstream cells flow through each cell. Archaeological ditches and natural drainage features both accumulate high flow values. D8 direction is computed via vectorized numpy slicing. Accumulation uses numba JIT (cached at module level) or pure Python fallback. """ logger.info(f" → Accumulation d'écoulement (flow accumulation)...") t0 = time.time() output = vis_dir / f"{basename}_flow_acc.tif" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np nan_mask = shared.nan_mask filled = shared.filled else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) rows, cols = dem_np.shape # Sink filling — priority-flood (O(n log n), NaN-aware) dem_filled = _priority_flood(filled, nan_mask) logger.info(f" ✓ Sink filling terminé ({time.time()-t0:.1f}s)") # D8 flow direction — vectorized via numpy slicing 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, 1, mode='constant', constant_values=np.nanmax(dem_filled[~np.isnan(dem_filled)]) + 10000) for d in range(8): nx = 1 + dx8[d] ny = 1 + dy8[d] neighbor_elev = padded[ny:ny + rows, nx:nx + cols] slope = (dem_filled - neighbor_elev) / (dist8[d] * resolution) slope[nan_mask] = -1 better = slope > max_slope flow_dir[better] = d max_slope[better] = slope[better] logger.info(f" ✓ Direction D8 terminée ({time.time()-t0:.1f}s)") # D8 accumulation — numba JIT (module-level cache) or pure Python fallback result = _d8_accumulate_numba(dem_filled, flow_dir, nan_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 logger.info(f" Accumulation D8 via Python (installez numba pour accélérer)") flat_dem = dem_filled[~nan_mask].flatten() valid_indices = np.where(~nan_mask.flatten())[0] sort_order = valid_indices[np.argsort(-flat_dem)] flow_acc = np.ones((rows, cols), dtype=np.float32) flow_acc[nan_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 nan_mask[nr, nc]: flow_acc[nr, nc] += flow_acc[r, c] logger.info(f" ✓ D8 accumulation terminé ({time.time()-t0:.1f}s)") # Log transform flow_result = np.log1p(flow_acc) flow_result[nan_mask] = np.nan _save_tif(output, flow_result, transform, crs) logger.info(f" ✓ Flow accumulation terminé ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Erreur flow accumulation: {e}", exc_info=True) return None # ============================================================ # Anomaly Mask — automatic threshold detection # ============================================================ def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, n_sigma=2.0): """Composite anomaly mask — automatic threshold detection (GPU if available). Reads pre-computed visualization layers (MSRM, SVF, Wavelet, Openness Neg, Roughness), normalizes each to z-scores, and combines them into a composite anomaly score. Pixels beyond `n_sigma` standard deviations of the local mean are flagged as suspicious. The output is a continuous score (0–1) where: - 0 = no anomaly (flat/natural terrain) - 1 = high anomaly (potential archaeological structure) This mask is directly usable in GIS for polygon extraction and field survey planning. Args: n_sigma: Number of standard deviations for the anomaly threshold. Lower = more sensitive (more false positives). Default 2.0 (good balance for archaeological detection). """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Détection automatique d'anomalies (seuil {n_sigma}σ){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_anomaly.tif" try: if shared: transform = shared.transform crs = shared.crs nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) rows, cols = nan_mask.shape # Collect available visualization layers from disk # Each is loaded, converted to |z-score|, and contributes to a weighted sum. # The weighted sum of |z| acts as a "vote": pixels where multiple layers # show anomalies get higher scores than pixels where only one layer fires. layer_configs = [ # (filename_pattern, weight) ("mslrm", 2.5), # Multi-scale relief — strongest signal ("negative_openness", 2.0), # Fossés, dolines ("roughness", 1.8), # Surface irregularity ("wavelet", 1.5), # Circular + linear structures ("svf", 1.3), # Sky-view depressions ("flow_acc", 1.2), # Drainage channels / ditches ("positive_openness", 0.8), # Surélevations ] layers = [] for pattern, weight in layer_configs: layer_path = vis_dir / f"{basename}_{pattern}.tif" if not layer_path.exists(): layer_path = vis_dir / f"{basename}_negative_openness.tif" if "neg" in pattern else None if layer_path is None or not layer_path.exists(): continue try: with rasterio.open(layer_path) as src: data = src.read(1).astype(np.float64) # Absolute z-score (captures both positive and negative deviations) valid = data[~nan_mask] if len(valid) == 0: continue mean_val = np.nanmean(valid) std_val = max(np.nanstd(valid), 0.01) zscore = np.abs(data - mean_val) / std_val zscore[nan_mask] = 0.0 layers.append((zscore, weight)) except Exception as e: logger.debug(f" Couche {pattern} non disponible: {e}") continue if not layers: # Fallback: use MSRM + roughness computed on the fly logger.info(" Aucune couche trouvée — calcul MSRM + rugosité en direct...") dem_np_safe = shared.dem_np if shared else dem_np # Quick MSRM (single scale 10m for speed) sigma_px = max(5, 10.0 / resolution) local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px) if shared else \ _filter_nanaware(dem_np_safe, xp_gaussian_filter, sigma=sigma_px) quick_mslrm = np.abs(dem_np_safe - local_mean) quick_mslrm[nan_mask] = np.nan valid = quick_mslrm[~nan_mask] std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 quick_mslrm = quick_mslrm / std_val layers.append((quick_mslrm, 2.5)) # Quick roughness fine_size = max(3, int(3 / resolution)) if fine_size % 2 == 0: fine_size += 1 if shared: fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size) fine_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=fine_size) else: fine_mean = _filter_nanaware(dem_np_safe.astype(np.float64), xp_uniform_filter, size=fine_size) fine_mean_sq = _filter_nanaware(dem_np_safe.astype(np.float64)**2, xp_uniform_filter, size=fine_size) roughness = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0)) roughness[nan_mask] = np.nan valid_r = roughness[~nan_mask] std_val_r = max(np.nanstd(valid_r), 0.01) if len(valid_r) > 0 else 0.01 roughness = roughness / std_val_r layers.append((roughness, 1.5)) # Weighted SUM of |z-scores| (not RMS — each layer votes independently) combined = np.zeros((rows, cols), dtype=np.float64) total_weight = 0.0 for layer_data, weight in layers: combined += layer_data * weight total_weight += weight if total_weight > 0: combined = combined / total_weight # Adaptive threshold: suppress pixels below the (100 - n_sigma*10)th percentile. # With n_sigma=2.0 → 80th percentile: keep only the top 20% of signal. # This adapts to each tile's terrain instead of a fixed z-score cutoff. threshold_pct = max(50, min(95, 100 - n_sigma * 10)) combined_valid = combined[~nan_mask] if combined_valid.size == 0: logger.warning(" ✗ Détection anomalies : tuile 100 % nodata") return None threshold_val = np.percentile(combined_valid, threshold_pct) combined = np.clip(combined - threshold_val, 0, None) # Rescale survivors to 0–1 above_thresh = combined[combined > 0] if len(above_thresh) > 0: p95 = np.percentile(above_thresh, 95) if p95 > 0: combined = np.clip(combined / p95, 0, 1) combined[nan_mask] = np.nan _save_tif(output, combined.astype(np.float32), transform, crs) n_anomaly = int(np.sum(combined > 0.1)) if np.any(combined > 0) else 0 pct_anomaly = n_anomaly / max(np.sum(~nan_mask), 1) * 100 logger.info(f" ✓ Détection anomalies terminée ({time.time()-t0:.1f}s) — " f"{pct_anomaly:.1f}% de la zone ({n_anomaly} px) au-delà de {n_sigma}σ") return output except Exception as e: logger.error(f" ✗ Erreur détection anomalies: {e}", exc_info=True) return None