"""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 15 m kernel (reused by 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(" → Computing filled DEM (NaN interpolation)...") 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(" → Computing 15 m LRM...") 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(" → Computing 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" GPU gradient failed ({e}) — falling back to 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. Uses a distance transform (O(n), vectorized): the indices of the nearest valid neighbor come out in a single pass. NearestNDInterpolator built a cKDTree over ALL valid points (25 M at 0.2 m) — several seconds per tile with holes, paid on the first access to SharedDEM.filled. """ nan_mask = np.isnan(arr) if not np.any(nan_mask): return arr, nan_mask # GPU (cupyx): the distance transform over 25 M pixels costs ~3 s on CPU, # the longest preparation step once everything else runs on 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" GPU hole filling failed ({e}) — falling back to 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. Optimization: the TANGENT of the angle (dz/dist) is accumulated instead of the angle itself — atan being strictly increasing, max(angles) = atan(max(tangents)). The (costly, full-image) arctan is therefore only applied at radius checkpoints, not at every ray step. A single pair of running maxima is kept and snapshotted at each checkpoint (radii are nested; previously each checkpoint redid the same computation). fmax ignores the padding NaNs: no more 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 sorted by step: (step, r_idx). Snapshots are taken when the # current step reaches the checkpoint step — rays therefore never go past # the largest requested checkpoint (equivalent to the former break). 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 (up to the last 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)) # Running max of the tangents (float32: half the VRAM of 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 # Angle tangents (positive: terrain above, negative: below). fmax # keeps the non-NaN operand: the padding border does not count, # as with the former where(isnan) — in a single operation. 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 of the reached checkpoint: converted to an angle ONCE 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 never reached (invalid steps): final state of the sweep 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 with automatic CPU fallback when VRAM runs out. Multi-radius 0.2 m tiles (5000×5000 px) can exceed the available VRAM (GPU shared with other services): rather than giving up on the visualization, the GPU is disabled for this worker and the computation is rerun on 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(" ⚠ Not enough VRAM (ray-tracing) — falling back to CPU for this 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" → Multi-directional hillshade{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 # No GPU copy of the raw DEM here: only gradient/slope/aspect are # used (already shared) — saves ~100 MB of VRAM per tile 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) # Spacing = resolution (m/px): without it the slope is wrong 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 done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Hillshade error: {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" → 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) # Spacing = resolution (m/px): without it the slope is wrong 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" ✓ Slope done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Slope error: {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 (slope orientation){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) # Spacing = resolution (m/px): without it the slope is wrong 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 done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Aspect error: {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 (multi-radius ray-tracing){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) # true radius in pixels, not truncated 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 done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ SVF error: {e}", exc_info=True) return None # Downsampling of the openness computation: ray-tracing is the most # expensive step of the pipeline (100 m radius = 500 steps at 0.2 m/px). # Openness being a smooth angular field (mean of horizons up to 100 m), it is # computed on a block-decimated grid — local max for positive openness, min # for negative, to preserve the relief that bounds the horizon angle (banks, # walls) — then resampled to the requested resolution. Cost ÷ factor³ (cells # ÷ factor², ray steps ÷ factor); factor 2 ≈ 8× faster, near-identical # rendering. 1 = full resolution. OPENNESS_DOWNSAMPLE = 2 # Ray-tracing radii (meters), averaged with equal weights. OPENNESS_RADII_M = (25, 50, 100) # Openness normalization references (degrees: mean, standard deviation). # A per-tile z-score made the scale non-seamless — same physical openness, # different color from one tile to the next depending on the surrounding # relief. Cross-tile medians measured at 0.2 m (×2 decimation) on 15 tiles # spread across the territory (plain, bocage, forest, mountain, volcanoes, # delta, urban, coast). FROZEN references: same openness = same color on # every tile. 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()} (multi-radius ray-tracing){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 # Block-decimated grid (see 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" Grid decimated ×{factor} ({rows}×{cols}) — resampled at the end") max_dist = int(max(radii_m) / res) # true radius in pixels, not truncated 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) # Back to the full-resolution grid (bilinear; the edge is padded by # edge replication if the dimensions were not divisible by the factor) 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 # Deviation from the frozen references, in sigmas: same physical # openness = same value on every tile → seamless color mosaic 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} done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Openness error: {e}", exc_info=True) return None # ============================================================ # Oriented relief: local openness × aspect in a single RGB image # ============================================================ # # Lightness (CIELAB L*) = micro-relief: mean positive horizon angle computed # on the detrended DTM (DTM − 10 m Gaussian), short radii 5/10/20 m (low # angle = open = light), plus a light directional shading. Hue = slope # orientation (aspect) on the CIELAB hue circle, at constant lightness: no # color creates false relief. No per-tile statistic (FIXED log scale) and a # total support (20 m radius + 4σ = 40 m detrend window, ~60 m in all) # narrower than the 100 m edge buffer: adjacent tiles join without seams. # # Cost kept under control: detrending and ray-tracing run on a grid # decimated to ~0.8 m (block max for the rays, block mean for the trend) — # only the final openness is upsampled; a dedicated kernel accumulates only # the mean of the angles (CuPy RawKernel on GPU, parallel numba on CPU, # numpy as a last resort); colors come from a precomputed table (L* × hue) # instead of a per-pixel Lab → sRGB conversion. RELIEF_DETREND_M = 10.0 # σ of the Gaussian trend removed from the DTM RELIEF_RADII_M = (5, 10, 20) # ray-tracing radii, equal weights RELIEF_N_DIRS = 16 # 16 directions: no 8-pointed star artifacts RELIEF_GRID_M = 0.8 # step of the decimated computation grid RELIEF_OPEN_RANGE = (0.5, 12.0) # degrees, fixed log scale (seamless) RELIEF_CHROMA = 60.0 # max CIELAB chroma (reached around L* = 50) RELIEF_SHADE_WEIGHT = 0.35 # share of directional shading in L* # Shading azimuth, expressed in the aspect frame used below # (arctan2(dy, dx), dy pointing south) — tuned on the rendering gallery. 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], clipped.""" 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 L* levels, 360 hues) → sRGB uint8, computed once. Chroma decreases toward black and white (L*(100−L*)/2500): colors stay inside the sRGB gamut instead of being clipped. """ 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): """(col, row) offsets of every ray step, distances and checkpoints. Same geometry as _ray_trace_horizons_core: direction k at angle 2πk/n, steps rounded to the pixel, distance = step × resolution. """ 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): """One CUDA thread per pixel; everything stays in registers (no intermediate array per direction or per radius).""" 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): """Same kernel in parallel numba (one pixel row per CPU thread).""" 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): """Vectorized fallback (no numba, no GPU): one array shift per step.""" 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): """Mean (directions × radii) of the positive horizon angle, in radians. GPU (CuPy RawKernel) if available — the result then stays on the GPU —, otherwise numba, otherwise numpy. A GPU failure (NVRTC compilation, VRAM) switches this computation to CPU without losing the visualization. Returns (angle, engine name). """ 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" ⚠ Oriented-relief GPU kernel unavailable ({e}) — falling back to 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): """Pad by edge replication (dimensions not a multiple of the factor).""" 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(): """Scalar rendering constants (shading without per-pixel trigonometry). cos(slope) = 1/√(1+g²) and sin(slope)·cos(az − aspect) = (cos az·dx + sin az·dy)/√(1+g²), with aspect = arctan2(dy, dx): shading costs only one square root per 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): """Vectorized colorization (CuPy on GPU, numpy as fallback).""" 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): """Single-pass CPU colorization (parallel numba): openness is read from the decimated grid by bilinear interpolation aligned on pixel centers (equivalent to xp_zoom), with no intermediate full-resolution array.""" 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): """Oriented relief: a single RGB image merging local openness and aspect. L* (lightness) = local positive openness (detrended DTM, short radii) 65 % + directional shading 35 %; hue = slope orientation; fixed CIELAB chroma. Output: RGB uint8 GeoTIFF (rendered as is, like the IGN base layers). """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Oriented relief (local openness × 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. Decimated grid: block max (preserves the relief that bounds the # horizon) and block mean (support of the trend to remove) 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. Detrending (DTM − Gaussian) on the decimated grid 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. Local openness: mean horizon angle (radians) on the decimated grid 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. Colors: openness upsampling, shading, L* × hue table 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" ✓ Oriented relief done ({time.time()-t0:.1f}s: " f"preparation {t_prep:.1f}s, rays {t_ray:.1f}s [{engine}], " f"colors {t_color:.1f}s [{engine_color}])") return output except Exception as e: logger.error(f" ✗ Oriented relief error: {e}", exc_info=True) return None # Ground point density: 16 levels on a fixed log scale (same gray = same # density everywhere, seamless mosaic). Level k starts at # DENSITY_BASE × 2^(k/2) pts/m²: 0 = < 0.35 (black, including no point at # all), 15 = ≥ 45 (white); two levels up = density doubled. DENSITY_LEVELS = 16 DENSITY_BASE = 0.25 def density_levels(density): """Level 0..15 of a density (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): """Density of the ground points kept for the DTM, in 16 gray levels. Read from the sidecar file written during rasterization (dtm.density_path): 1 m grid, kept as is (density carries no finer detail — an image 25× lighter than at 0.2 m). Output: level 0..15 (float32). """ logger.info(" → Ground point density...") 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" ✗ Density file missing ({src_path.name}): regenerate the 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" ✓ Point density done ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Point density error: {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 scales: {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 done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ MSRM error: {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 (adaptive LRM){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 # Per-tile z-score: units comparable between tiles, fixed ±3σ # rendering range → homogeneous color mosaic 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 done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ SAILORE error: {e}", exc_info=True) return None # ============================================================ # Roughness # ============================================================ def _integral_sums(x): """2D integral sums: S[i,j] = sum of x[0:i, 0:j] (float64). Row/column 0 are zeros — gives the sum of any window from 4 corners, including against the border (index 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): """Local standard deviation over a size×size window via integral sums. Cost independent of the window size (4 lookups per pixel) — versus a uniform_filter whose cost grows with the window (75 px at 0.2 m for the broad scale). At the borders the window is truncated and normalized by the actual number of elements (tiles overlap through the edge buffer, so the border has no visible effect). """ 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) # Actual number of elements in the window (truncated at the borders) 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)) # Roughness normalization references (meters of local standard deviation). # Cross-tile medians measured on 20 real tiles at 0.2 m: per-tile # normalization (z-score) made the scale non-seamless — the deviation ranged # from 0.05 to 0.58 m depending on the tile for the fine scale. FROZEN # references: same physical roughness = same value on every tile. 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. Local standard deviations are computed from integral sums: two cumsums shared by both scales, 4-corner extraction — constant cost whatever the window (a uniform_filter costs in proportion to its size, 75 px at 0.2 m for the broad scale). Normalized by FIXED references (measured medians), not per tile: mosaics are seamless, same value = same color. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Surface roughness{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) # Integral sums shared by both scales (X and 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 # Weighted combination on a common physical scale (fixed references: # seamless between tiles — see module constants) 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" ✓ Roughness done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Roughness error: {e}", exc_info=True) return None # ============================================================ # Surface exposure (solar illumination) # ============================================================ 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(" → Surface exposure (solar illumination)...") 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" ✓ Surface exposure done ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Surface exposure error: {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 (scales below max(2 × resolution, 1 m) are dropped, so the 0.5 m candidate is never used). 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" → Multi-scale Mexican Hat wavelet{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) # 100 m dropped: it mostly responds to large landforms, not to # structures. Emphasis on 1-10 m (small 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" CWT scales: {scales} m (resolution {resolution} m/px)") from scipy.ndimage import gaussian_laplace, gaussian_filter # Removal of large landforms (hills, valleys): a Gaussian local mean # is subtracted before the CWT. Gaussian smoothing preserves planar # slopes, so the residual is analyzed relative to its ~35 m # neighborhood: a ditch on a hilltop or hillside does not stand out # more than the same ditch on flat ground. # Fraction kept for a Gaussian structure of width σ_f: # σ_t²/(σ_f²+σ_t²) → 10 m: 92%, 20 m: 75%, 150 m hill: 5%. # Measured on a noisy synthetic DTM: the contrast of small structures # is insensitive to σ_t; only the hilltop/flat background ratio # changes (1.71 without detrending → 1.21 at 35 m). 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" Large landforms removed (Gaussian trend {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] # Robust σ (MAD): the classic std is inflated by the tails # (strong structures, tile edges) and varies a lot from one tile # to the next — the main cause of per-tile color casts on the # map. 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 # Rescaling by the tile median: the RMS becomes an index relative to # the tile's typical level (median = 1). The distribution is then # comparable from one tile to the next — a prerequisite for a single # fixed color stretch that is homogeneous across tiles. 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" ✓ Wavelet done ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}") return output except Exception as e: logger.error(f" ✗ Wavelet error: {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. Returns None if numba is unavailable (the caller falls back to Python). """ 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(" → 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 done ({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" ✓ D8 direction done ({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(" D8 accumulation via numba") else: # Pure Python fallback logger.info(" D8 accumulation via Python (install numba to speed it up)") 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 done ({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 done ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Flow accumulation error: {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, Flow accumulation, Openness Pos), normalizes each to absolute z-scores, and combines them into a weighted composite anomaly score. Pixels below the (100 − 10 × n_sigma)th percentile of the score (clamped to 50–95) are zeroed; the rest 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: Threshold control (percentile = 100 − 10 × n_sigma). Lower = more sensitive (more false positives). Default 2.0 (80th percentile, top 20% of the score kept). """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → Automatic anomaly detection (threshold {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), # Ditches, sinkholes ("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), # Raised features ] 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" Layer {pattern} unavailable: {e}") continue if not layers: # Fallback: use MSRM + roughness computed on the fly logger.info(" No layer found — computing MSRM + roughness on the fly...") 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(" ✗ Anomaly detection: tile is 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" ✓ Anomaly detection done ({time.time()-t0:.1f}s) — " f"{pct_anomaly:.1f}% of the area ({n_anomaly} px) above the {n_sigma}σ threshold") return output except Exception as e: logger.error(f" ✗ Anomaly detection error: {e}", exc_info=True) return None