"""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 time import warnings from pathlib import Path import numpy as np import rasterio from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, xp_minimum_filter, xp_maximum_filter, gpu_cleanup from . import gpu as _gpu_mod logger = logging.getLogger("lidar") # CuPy module reference — lazily imported on first GPU use. # If disable_gpu() is called at runtime, HAS_GPU becomes False # and xp delegates to numpy instead. _cp = None class _XPProxy: """Proxy that delegates array operations to cupy or numpy. Checks HAS_GPU on every attribute access so that disable_gpu() (called on CUDA errors) takes effect immediately, without needing to change every call site in visualizations.py. """ def __getattr__(self, name): global _cp from . import gpu as _gpu_mod if _gpu_mod.HAS_GPU: if _cp is None: try: import cupy _cp = cupy except ImportError: pass if _cp is not None: return getattr(_cp, name) return getattr(np, name) xp = _XPProxy() class SharedDEM: """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, aspect, curvature) - LRM at 15m kernel (shared by lrm + anomalies) 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, aspect) 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 aspect(self): self._ensure_gradient() return self._gradient[4] @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) aspect = np.mod(np.degrees(np.arctan2(dy, dx)), 360) self._gradient = (dy, dx, slope_rad, slope_deg, aspect) @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. """ from scipy.interpolate import NearestNDInterpolator nan_mask = np.isnan(arr) if not np.any(nan_mask): return arr, nan_mask valid = ~nan_mask y_coords, x_coords = np.where(valid) if len(y_coords) == 0: return np.zeros_like(arr), nan_mask z_values = arr[valid] interp = NearestNDInterpolator( np.column_stack((y_coords, x_coords)), z_values ) y_missing, x_missing = np.where(nan_mask) filled = arr.copy() filled[y_missing, x_missing] = interp(y_missing, x_missing) 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 # ============================================================ # 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" 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_tag}") 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_tag}") 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.""" gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Aspect (Orientation){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_aspect.tif" try: if shared: transform = shared.transform crs = shared.crs aspect = shared.aspect nan_mask = shared.nan_mask if _gpu_mod.HAS_GPU: aspect = to_gpu(aspect) else: dem_np, transform, crs = _read_dem(dem_file) dem = to_gpu(dem_np) dy, dx = xp.gradient(dem) aspect = xp.arctan2(dy, dx) * 180 / xp.pi aspect = xp.mod(aspect, 360) nan_mask = np.isnan(dem_np) _save_tif(output, to_cpu(aspect) if _gpu_mod.HAS_GPU else aspect, transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur aspect: {e}", exc_info=True) return None def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): """Generate curvature (terrain concavity/convexity) map — GPU if available.""" gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Courbure (Curvature){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_curvature.tif" try: if shared: transform = shared.transform crs = shared.crs dx = shared.dx dy = shared.dy nan_mask = shared.nan_mask if _gpu_mod.HAS_GPU: dx = to_gpu(dx) dy = to_gpu(dy) 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) d2z_dx2 = xp.gradient(dx, axis=1) d2z_dy2 = xp.gradient(dy, axis=0) curvature = (d2z_dx2 + d2z_dy2) / 2 _save_tif(output, to_cpu(curvature), transform, crs, nan_mask=nan_mask) logger.info(f" ✓ Courbure terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur curvature: {e}", exc_info=True) return None # ============================================================ # GPU-accelerated visualizations # ============================================================ def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None): """Local Relief Model - deviation from local mean (GPU if available). Kernel sigma adapts to resolution: finer kernel at higher resolution to capture micro-relief details. At 0.5m/px: 15m, at 0.2m/px: ~5m. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Local Relief Model{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_lrm.tif" try: if shared: transform = shared.transform crs = shared.crs lrm = shared.lrm_15.copy() else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) # Adapt sigma to resolution: standard 15m at 0.5m, finer at higher res sigma_m = max(5.0, 15.0 * 0.5 / resolution) logger.info(f" LRM sigma={sigma_m:.1f}m (résolution {resolution}m/px)") local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_m / resolution) lrm = dem_np - local_mean lrm[nan_mask] = np.nan _save_tif(output, lrm.astype(np.float32), transform, crs) logger.info(f" ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur LRM: {e}", exc_info=True) return None def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): """Sky-View Factor - ray-tracing on 16 azimuths (GPU if available). For each pixel, trace rays in N directions, find the max horizon angle in each direction, then SVF = (1/N) * sum(cos²(horizon_angle)). Valleys/crevices have low SVF (obstructed sky), ridges/peaks have high SVF. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Sky-View Factor (ray-tracing){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_svf.tif" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) rows, cols = dem_np.shape res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 16 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx_dir = np.cos(angles) dy_dir = np.sin(angles) # Cap max_dist to avoid excessive computation at high resolution # 100m radius is sufficient; at 0.2m that's 500 steps which is very slow max_dist = min(int(100 / res), 300) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) svf = xp.zeros_like(dem) for d_idx in range(n_dirs): ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] horizon = xp.zeros_like(dem) # Pre-compute all valid steps for this direction valid_steps = [] for step in range(1, max_dist + 1): px = int(round(ddx * step)) py = int(round(ddy * step)) dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) if dist_m < res * 0.5: continue valid_steps.append((step, px, py, dist_m)) # Batch all shifts into a single array for vectorized max computation for step, px, py, dist_m in valid_steps: elev_diff = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] - dem angle = xp.arctan2(elev_diff, dist_m) horizon = xp.where(xp.isnan(angle), horizon, xp.maximum(horizon, xp.nan_to_num(angle, nan=0))) # SVF uses cos²(horizon angle) — fraction of visible sky svf += xp.cos(horizon) ** 2 svf /= n_dirs svf_np = to_cpu(svf).astype(np.float32) svf_np[nan_mask] = np.nan _save_tif(output, svf_np, transform, crs) logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur SVF: {e}", exc_info=True) return None def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None): """Positive/Negative Openness - true zenith/nadir angle computation (GPU if available). For each pixel, in 8 directions (N, NE, E, SE, S, SW, W, NW): - Positive openness: max zenith angle (angle from vertical to highest visible terrain) - Negative openness: max nadir angle (angle from vertical down to lowest terrain) Result is averaged across all 8 directions. Ray radius adapts to resolution: 100m for better detection of large enclosures. """ 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){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_{name}.tif" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) rows, cols = dem_np.shape res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx_dir = np.cos(angles) dy_dir = np.sin(angles) max_dist = min(int(100 / res), 300) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) openness_sum = xp.zeros_like(dem) for d_idx in range(n_dirs): ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] max_angle = xp.zeros_like(dem) for step in range(1, max_dist + 1): px = int(round(ddx * step)) py = int(round(ddy * step)) dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) if dist_m < res * 0.5: continue elev_diff = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] - dem if positive: angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) else: angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) max_angle = xp.where(xp.isnan(angle), max_angle, xp.maximum(max_angle, xp.nan_to_num(angle, nan=0))) openness_sum += max_angle openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32) openness_result[nan_mask] = np.nan _save_tif(output, openness_result, transform, crs) logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur openness: {e}", exc_info=True) return None def generate_mslrm(dem_file, basename, vis_dir, resolution, 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) candidate_scales = [2, 5, 10, 20, 50, 100, 200] sigmas = [s for s in candidate_scales if s >= min_scale] # Archaeological weights: favor 5-25m range (ditches, enclosures, tumulus) scale_weights = { 2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6, 200: 0.4, } 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: x / std — preserves sign and contrast better than z-score valid_lrm = lrm[~nan_mask] lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 lrm = lrm / lrm_std 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_tag}") return output except Exception as e: logger.error(f" ✗ Erreur MSRM: {e}", exc_info=True) return None def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None): """Multi-Scale Topographic Position Index (GPU if available). TPI = elevation - mean(neighborhood). Computed at 4 scales with std normalization and weighted combination. Weights favor fine and medium scales (archaeologically relevant). """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → TPI multi-échelle{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_tpi.tif" 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) # 4 scales: fine (3m), medium (15m), broad (50m), landscape (200m) scales_m = [3, 15, 50, 200] weights = [1.5, 2.0, 1.2, 0.5] # Favor medium scales (ditches, enclosures) tpi_stack = [] for scale_m, weight in zip(scales_m, weights): size = max(3, int(scale_m / resolution)) if size % 2 == 0: size += 1 if shared: local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=size) else: local_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=size) tpi = dem_np - local_mean tpi[nan_mask] = np.nan # Std normalization — preserves sign and contrast better than z-score valid = tpi[~nan_mask] tpi_std = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 tpi = tpi / tpi_std tpi_stack.append(tpi.astype(np.float32)) # Weighted combination tpi_array = np.array(tpi_stack) weights_3d = np.array(weights)[:, np.newaxis, np.newaxis] with np.errstate(invalid='ignore', divide='ignore'): with warnings.catch_warnings(): warnings.filterwarnings('ignore', message='Mean of empty slice') tpi_combined = np.nansum(tpi_array * weights_3d, axis=0) / np.sum(weights) tpi_combined[nan_mask] = np.nan _save_tif(output, tpi_combined.astype(np.float32), transform, crs) logger.info(f" ✓ TPI terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur TPI: {e}", exc_info=True) return None # ============================================================ # 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. """ 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 sigma_mid = (sigma_min + sigma_max) / 2 slope_norm = np.clip(slope_deg / 30.0, 0, 1) 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 if shared: lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) else: lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) lrm_medium[nan_mask] = np.nan 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 _save_tif(output, sailore.astype(np.float32), transform, crs) logger.info(f" ✓ SAILORE terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur SAILORE: {e}", exc_info=True) return None # ============================================================ # Roughness # ============================================================ def generate_roughness(dem_file, basename, vis_dir, resolution, 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. """ 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 else: dem_np, transform, crs = _read_dem(dem_file) nan_mask = np.isnan(dem_np) # Fine roughness (3m window) 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) fine_mean_sq[shared.nan_mask] = np.nan else: fine_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=fine_size) fine_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=fine_size) roughness_fine = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0)) roughness_fine[nan_mask] = np.nan # Broad roughness (15m window) broad_size = max(3, int(15 / resolution)) if broad_size % 2 == 0: broad_size += 1 if shared: broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size) broad_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=broad_size) broad_mean_sq[shared.nan_mask] = np.nan else: broad_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=broad_size) broad_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=broad_size) roughness_broad = np.sqrt(np.maximum(broad_mean_sq - broad_mean * broad_mean, 0)) roughness_broad[nan_mask] = np.nan # Std normalization per scale then weighted combination fine_valid = roughness_fine[~nan_mask] broad_valid = roughness_broad[~nan_mask] fine_std = max(np.nanstd(fine_valid), 0.01) if len(fine_valid) > 0 else 0.01 broad_std = max(np.nanstd(broad_valid), 0.01) if len(broad_valid) > 0 else 0.01 roughness = 0.7 * roughness_fine / fine_std + 0.3 * roughness_broad / broad_std roughness[nan_mask] = np.nan _save_tif(output, roughness, transform, crs) logger.info(f" ✓ Rugosité terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur rugosité: {e}", exc_info=True) return None # ============================================================ # Wavelet # ============================================================ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): """Mexican Hat wavelet multi-scale analysis (GPU if available). CWT 2D at multiple scales adapted to resolution. - At 0.5m/px: [1, 2, 5, 10, 20, 50, 100]m - At 0.2m/px: [0.5, 1, 2, 5, 10, 20, 50, 100]m - Higher resolution = more fine scales available Uses std normalization per scale and weighted combination with emphasis on archaeologically relevant scales (2-50m). """ 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)) # Adapt scales to resolution: finer scales available at higher resolution min_scale = max(resolution * 2, 1.0) candidate_scales = [0.5, 1, 2, 5, 10, 20, 50, 100] scales = [s for s in candidate_scales if s >= min_scale] # Weights favor archaeological scales (2-50m: ditches, enclosures, tumulus) scale_weights = { 0.5: 0.6, # Fine texture 1.0: 0.8, # Micro-relief 2.0: 1.5, # Small ditches, paths — key scale 5.0: 2.0, # Fossés, small enclosures — key archaeological scale 10.0: 1.8, # Medium structures 20.0: 1.5, # Large enclosures, tumulus 50.0: 1.0, # Very large enclosures 100.0: 0.6, # Landscape-level features } 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)") 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(filled), sigma=sigma_px) response = to_cpu(response) except Exception: from scipy.ndimage import gaussian_laplace response = -gaussian_laplace(filled, sigma=sigma_px) else: from scipy.ndimage import gaussian_laplace response = -gaussian_laplace(filled, sigma=sigma_px) response[nan_mask] = np.nan # Std normalization: scale by standard deviation to make scales comparable valid = response[~nan_mask] std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 response = response / std_val wavelet_stack.append(response) # Weighted RMS: sqrt(sum(w * x²) / sum(w)) # Preserves contrast at key archaeological scales 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 _save_tif(output, combined.astype(np.float32), transform, crs) logger.info(f" ✓ Ondelette terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur ondelette: {e}", exc_info=True) return None # ============================================================ # Anisotropic Openness # ============================================================ # Path Detection (chemins et sentiers) # ============================================================ def generate_paths(dem_file, basename, vis_dir, resolution, shared=None): """Cheminement — openness directionnelle maximale pour détecter chemins et sentiers. Pour chaque direction (8 directions), calcule openness positive - négative, puis prend le maximum sur toutes les directions. Les chemins et sentiers ressortent en valeurs élevées quelle que soit leur orientation. Contrairement à l'openness anisotropique qui privilégie NW-SE et NE-SW, cette visualisation traite toutes les directions de manière égale et combine positive et négative en une seule image. Les chemins perpendiculaires à une direction auront une forte différence dans cette direction, donc le maximum sur toutes les directions les fait ressortir. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Cheminement (chemins et sentiers){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_paths.tif" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) rows, cols = dem_np.shape res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx_dir = np.cos(angles) dy_dir = np.sin(angles) max_dist = min(int(100 / res), 300) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) max_diff = xp.full_like(dem, -1e6) # Will track max over all directions for d_idx in range(n_dirs): ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] # Positive openness: max zenith angle in this direction max_pos_angle = xp.zeros_like(dem) # Negative openness: max nadir angle in this direction max_neg_angle = xp.zeros_like(dem) for step in range(1, max_dist + 1): px = int(round(ddx * step)) py = int(round(ddy * step)) dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) if dist_m < res * 0.5: continue elev_diff = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] - dem # Positive: angle to terrain above viewer pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) max_pos_angle = xp.where(xp.isnan(pos_angle), max_pos_angle, xp.maximum(max_pos_angle, xp.nan_to_num(pos_angle, nan=0))) # Negative: angle to terrain below viewer neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) max_neg_angle = xp.where(xp.isnan(neg_angle), max_neg_angle, xp.maximum(max_neg_angle, xp.nan_to_num(neg_angle, nan=0))) # Difference highlights linear features perpendicular to this direction diff = max_pos_angle - max_neg_angle max_diff = xp.maximum(max_diff, diff) paths_result = to_cpu(max_diff).astype(np.float32) paths_result[nan_mask] = np.nan _save_tif(output, paths_result, transform, crs) logger.info(f" ✓ Cheminement terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur cheminement: {e}", exc_info=True) return None # ============================================================ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None): """Anisotropic Openness - weighted directional openness emphasizing oblique directions (GPU if available). Computes positive and negative openness with anisotropic weighting: NW/SE directions weighted more heavily to enhance detection of structures aligned NE-SW (common in French archaeological sites: villas, enclosures). The anisotropic weighting makes subtle linear features more visible than standard isotropic openness which averages all directions equally. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Openness Anisotropique{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_aniso_open.tif" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np rows, cols = dem_np.shape res = resolution dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask else: dem_np, transform, crs = _read_dem(dem_file) rows, cols = dem_np.shape res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx_dir = np.cos(angles) dy_dir = np.sin(angles) # Anisotropic weights: emphasize NW-SE and NE-SW directions # These orientations are most productive for detecting archaeological features # aligned with Roman and medieval settlement patterns in France weights = np.array([1.0, 1.5, 1.0, 1.5, 1.0, 1.5, 1.0, 1.5]) max_dist = min(int(100 / res), 300) padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) pos_sum = xp.zeros_like(dem) neg_sum = xp.zeros_like(dem) weight_total = 0.0 for d_idx in range(n_dirs): ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] w = weights[d_idx] weight_total += w max_pos_angle = xp.zeros_like(dem) max_neg_angle = xp.zeros_like(dem) for step in range(1, max_dist + 1): px = int(round(ddx * step)) py = int(round(ddy * step)) dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) if dist_m < res * 0.5: continue elev_diff = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] - dem # Positive openness: max zenith angle pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) max_pos_angle = xp.where(xp.isnan(pos_angle), max_pos_angle, xp.maximum(max_pos_angle, xp.nan_to_num(pos_angle, nan=0))) # Negative openness: max nadir angle neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) max_neg_angle = xp.where(xp.isnan(neg_angle), max_neg_angle, xp.maximum(max_neg_angle, xp.nan_to_num(neg_angle, nan=0))) pos_sum += max_pos_angle * w neg_sum += max_neg_angle * w # Combined: positive minus negative openness (anisotropic) pos_avg = pos_sum / weight_total neg_avg = neg_sum / weight_total aniso_result = to_cpu(xp.degrees(pos_avg - neg_avg)).astype(np.float32) aniso_result[nan_mask] = np.nan _save_tif(output, aniso_result, transform, crs) logger.info(f" ✓ Openness anisotropique terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur openness anisotropique: {e}", exc_info=True) return None