"""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 = 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 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 dem = to_gpu(shared.dem_np) 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) dy, dx = xp.gradient(dem) 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) 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 _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) dy, dx = xp.gradient(dem) 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 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). Results are combined with equal weight across radii, then normalized by standard deviation for cross-tile comparability. """ 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 = [25, 50, 100] max_dist = int(max(radii_m) / res) # rayon réel en pixels, non tronqué n_dirs = 8 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) openness_result[nan_mask] = np.nan # Z-score (écarts locaux en sigmas) : unités comparables entre tuiles, # plage de rendu fixe → mosaïque de couleur homogène valid = openness_result[~nan_mask] if len(valid) > 0: std_val = max(np.nanstd(valid), 0.01) openness_result = (openness_result - np.nanmean(valid)) / std_val _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 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)) threshold_val = np.percentile(combined[~nan_mask], 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