6780f2e25b
Remplace NearestNDInterpolator (cKDTree sur les 25 M de points valides, plusieurs secondes par dalle trouée en mode IGN pur) par distance_transform_edt : plus proche voisin en une passe O(n). Voisins équidistants arbitres differemment, sans effet sur les rendus.
1562 lines
62 KiB
Python
1562 lines
62 KiB
Python
"""Terrain visualization functions for LiDAR archaeological analysis.
|
||
|
||
Each function takes (dem_file, basename, vis_dir, resolution) as explicit
|
||
parameters and returns the path to the output GeoTIFF file, or None on error.
|
||
|
||
When a SharedDEM object is provided via the `shared` parameter, pre-computed
|
||
data (gradient, NaN mask, LRM) is reused across visualizations to avoid
|
||
redundant I/O and computation.
|
||
"""
|
||
|
||
import logging
|
||
import math
|
||
import time
|
||
import warnings
|
||
from pathlib import Path
|
||
|
||
import numpy as np
|
||
import rasterio
|
||
from .gpu import to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, gpu_cleanup
|
||
from . import gpu as _gpu_mod
|
||
|
||
logger = logging.getLogger("lidar")
|
||
|
||
# CuPy module reference — lazily imported on first GPU use.
|
||
# If disable_gpu() is called at runtime, HAS_GPU becomes False
|
||
# and xp delegates to numpy instead.
|
||
_cp = None
|
||
|
||
|
||
class _XPProxy:
|
||
"""Proxy that delegates array operations to cupy or numpy.
|
||
|
||
Checks HAS_GPU on every attribute access so that disable_gpu()
|
||
(called on CUDA errors) takes effect immediately, without needing
|
||
to change every call site in visualizations.py.
|
||
"""
|
||
|
||
def __getattr__(self, name):
|
||
global _cp
|
||
if _gpu_mod.HAS_GPU:
|
||
if _cp is None:
|
||
try:
|
||
import cupy
|
||
_cp = cupy
|
||
except ImportError:
|
||
pass
|
||
if _cp is not None:
|
||
return getattr(_cp, name)
|
||
return getattr(np, name)
|
||
|
||
|
||
xp = _XPProxy()
|
||
|
||
|
||
class SharedDEM:
|
||
"""Pre-computed DEM data shared across all visualizations.
|
||
|
||
Reads the DEM once and lazily computes on first access:
|
||
- NaN mask and filled DEM (avoids 20+ calls to _fill_nans)
|
||
- Gradient components (shared by hillshade, slope)
|
||
- LRM at 15m kernel (shared by mslrm + sailore)
|
||
|
||
Attributes are computed lazily on first access to avoid computing
|
||
data that is never used (e.g. LRM when only hillshade needs generation).
|
||
"""
|
||
|
||
def __init__(self, dem_file, resolution):
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
self.dem_file = dem_file
|
||
self.resolution = resolution
|
||
self.transform = transform
|
||
self.crs = crs
|
||
self.nan_mask = np.isnan(dem_np)
|
||
self.dem_np = dem_np.astype(np.float32)
|
||
|
||
# Lazy caches — computed on first access
|
||
self._filled = None
|
||
self._gradient = None # (dy, dx, slope_rad, slope_deg)
|
||
self._lrm_15 = None
|
||
|
||
# GPU lazy caches
|
||
self._filled_gpu = None
|
||
self._dem_gpu = None
|
||
|
||
@property
|
||
def filled(self):
|
||
"""Filled DEM (NaN interpolated) — computed lazily."""
|
||
if self._filled is None:
|
||
logger.debug(" → Calcul filled DEM (interpolation NaN)...")
|
||
self._filled, _ = _fill_nans(self.dem_np)
|
||
return self._filled
|
||
|
||
@property
|
||
def dy(self):
|
||
self._ensure_gradient()
|
||
return self._gradient[0]
|
||
|
||
@property
|
||
def dx(self):
|
||
self._ensure_gradient()
|
||
return self._gradient[1]
|
||
|
||
@property
|
||
def slope_rad(self):
|
||
self._ensure_gradient()
|
||
return self._gradient[2]
|
||
|
||
@property
|
||
def slope_deg(self):
|
||
self._ensure_gradient()
|
||
return self._gradient[3]
|
||
|
||
@property
|
||
def lrm_15(self):
|
||
"""LRM at 15m kernel — computed lazily."""
|
||
if self._lrm_15 is None:
|
||
logger.debug(" → Calcul LRM 15m...")
|
||
sigma_15 = 15.0 / self.resolution
|
||
local_mean_15 = _filter_nanaware_from_filled(self, xp_gaussian_filter, sigma=sigma_15)
|
||
self._lrm_15 = self.dem_np - local_mean_15
|
||
self._lrm_15[self.nan_mask] = np.nan
|
||
return self._lrm_15
|
||
|
||
def _ensure_gradient(self):
|
||
"""Compute gradient components lazily on first access."""
|
||
if self._gradient is None:
|
||
logger.debug(" → Calcul gradient...")
|
||
dy = np.gradient(self.filled, self.resolution, axis=0)
|
||
dx = np.gradient(self.filled, self.resolution, axis=1)
|
||
slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
|
||
slope_deg = np.degrees(slope_rad)
|
||
self._gradient = (dy, dx, slope_rad, slope_deg)
|
||
|
||
@property
|
||
def filled_gpu(self):
|
||
"""Lazy GPU copy of the filled DEM."""
|
||
if self._filled_gpu is None and _gpu_mod.HAS_GPU:
|
||
self._filled_gpu = to_gpu(self.filled)
|
||
return self._filled_gpu
|
||
|
||
@property
|
||
def dem_gpu(self):
|
||
"""Lazy GPU copy of the DEM."""
|
||
if self._dem_gpu is None and _gpu_mod.HAS_GPU:
|
||
self._dem_gpu = to_gpu(self.dem_np)
|
||
return self._dem_gpu
|
||
|
||
|
||
def _filter_nanaware_from_filled(shared, filter_func, *args, **kwargs):
|
||
"""Apply filter on pre-filled DEM data (skips expensive _fill_nans).
|
||
|
||
Uses the SharedDEM.filled array directly, then restores NaN mask.
|
||
If GPU is available, reuses the lazy GPU copy to avoid redundant transfers.
|
||
"""
|
||
if _gpu_mod.HAS_GPU:
|
||
filled_gpu = shared.filled_gpu
|
||
else:
|
||
filled_gpu = None
|
||
|
||
if filled_gpu is not None:
|
||
result_gpu = filter_func(filled_gpu, *args, **kwargs)
|
||
result = to_cpu(result_gpu)
|
||
gpu_cleanup()
|
||
else:
|
||
result = filter_func(shared.filled, *args, **kwargs)
|
||
result[shared.nan_mask] = np.nan
|
||
return result
|
||
|
||
|
||
def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None, nan_mask=None):
|
||
"""Helper to save a 2D or 3D array as GeoTIFF.
|
||
|
||
Args:
|
||
nan_mask: Optional boolean mask (True=NaN) to apply before saving.
|
||
Restores NaN zones in gradient-derived products that were
|
||
computed on the filled DEM.
|
||
"""
|
||
if nan_mask is not None:
|
||
data = np.array(data, dtype=dtype, copy=True)
|
||
data[nan_mask] = np.nan
|
||
|
||
# Auto-detect nodata for float types with NaN
|
||
if nodata is None and dtype.startswith('float') and np.any(np.isnan(data)):
|
||
nodata = float('nan')
|
||
|
||
if data.ndim == 2:
|
||
height, width = data.shape
|
||
with rasterio.open(
|
||
output_path, 'w', driver='GTiff',
|
||
height=height, width=width, count=count,
|
||
dtype=dtype, crs=crs, transform=transform,
|
||
compress='deflate', nodata=nodata
|
||
) as dst:
|
||
dst.write(data.astype(dtype), 1)
|
||
elif data.ndim == 3:
|
||
bands, height, width = data.shape
|
||
with rasterio.open(
|
||
output_path, 'w', driver='GTiff',
|
||
height=height, width=width, count=bands,
|
||
dtype=dtype, crs=crs, transform=transform,
|
||
compress='deflate', nodata=nodata
|
||
) as dst:
|
||
for i in range(bands):
|
||
dst.write(data[i].astype(dtype), i + 1)
|
||
|
||
|
||
def _read_dem(dem_file):
|
||
"""Read DEM file and return (data, transform, crs)."""
|
||
with rasterio.open(dem_file) as src:
|
||
return src.read(1), src.transform, src.crs
|
||
|
||
|
||
def _fill_nans(arr):
|
||
"""Fill NaN values using nearest-neighbor interpolation.
|
||
|
||
Returns (filled_array, nan_mask) so the caller can restore NaN after filtering.
|
||
|
||
Via transformée de distance (O(n), vectorisé) : les indices du plus proche
|
||
voisin valides sortent en une passe. NearestNDInterpolator construisait un
|
||
cKDTree sur TOUS les points valides (25 M à 0,2 m) — plusieurs secondes
|
||
par dalle trouée, payées au premier accès de SharedDEM.filled.
|
||
"""
|
||
nan_mask = np.isnan(arr)
|
||
if not np.any(nan_mask):
|
||
return arr, nan_mask
|
||
from scipy.ndimage import distance_transform_edt
|
||
_, (iy, ix) = distance_transform_edt(nan_mask, return_indices=True)
|
||
filled = arr.copy()
|
||
filled[nan_mask] = arr[iy[nan_mask], ix[nan_mask]]
|
||
return filled, nan_mask
|
||
|
||
|
||
def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
|
||
"""Apply a filter to an array while preserving NaN zones.
|
||
|
||
1. Fill NaN with nearest-neighbor interpolation
|
||
2. Apply the filter
|
||
3. Restore original NaN mask on the result
|
||
|
||
Args:
|
||
arr: Input array (numpy or cupy).
|
||
filter_func: Function that takes (array, *args, **kwargs) and returns filtered array.
|
||
use_gpu: If True, apply filter on GPU (send filled array to GPU first).
|
||
Returns:
|
||
Filtered array with original NaN positions preserved.
|
||
"""
|
||
is_gpu_arr = _gpu_mod.HAS_GPU and _cp is not None and isinstance(arr, _cp.ndarray)
|
||
arr_np = to_cpu(arr) if is_gpu_arr else arr
|
||
|
||
filled, nan_mask = _fill_nans(arr_np)
|
||
|
||
if use_gpu and _gpu_mod.HAS_GPU:
|
||
filled_gpu = to_gpu(filled)
|
||
result_gpu = filter_func(filled_gpu, *args, **kwargs)
|
||
result = to_cpu(result_gpu)
|
||
gpu_cleanup()
|
||
else:
|
||
result = filter_func(filled, *args, **kwargs)
|
||
|
||
result[nan_mask] = np.nan
|
||
return result
|
||
|
||
|
||
# ============================================================
|
||
# Shared ray-tracing core
|
||
# ============================================================
|
||
|
||
def _prepare_dem_for_raycast(dem_file, shared, resolution):
|
||
"""Load DEM and prepare padded array for ray-tracing.
|
||
|
||
Returns (dem_filled, dem_np, rows, cols, res, nan_mask,
|
||
transform, crs) ready for ray-tracing.
|
||
dem_filled is a CPU numpy array (filled, no NaN).
|
||
"""
|
||
if shared:
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem = shared.filled
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
filled, _ = _fill_nans(dem_np)
|
||
dem = filled
|
||
res = resolution
|
||
rows, cols = dem_np.shape
|
||
return dem, dem_np, rows, cols, res, nan_mask, transform, crs
|
||
|
||
|
||
def _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m=None):
|
||
"""Core ray-tracing: compute max zenith/nadir angles per direction and radius.
|
||
|
||
For each pixel, in each direction, traces rays outward up to max_dist steps,
|
||
recording the max upward angle (positive openness) and max downward angle
|
||
(negative openness) reached at each radius checkpoint.
|
||
|
||
Optimisation : on accumule la TANGENTE de l'angle (dz/dist) au lieu de
|
||
l'angle lui-même — atan étant strictement croissante, max(angles) =
|
||
atan(max(tangentes)). L'arctan (coûteuse, pleine image) n'est donc plus
|
||
appliquée qu'aux checkpoints de rayon, pas à chaque pas de rayon.
|
||
Un seul couple de max cumulés est maintenu, snapshoté à chaque checkpoint
|
||
(les rayons étant emboîtés, chaque checkpoint réutilisait avant le même
|
||
calcul 3 fois). fmax ignore les NaN du padding : plus de nan_to_num/where.
|
||
|
||
Padding on CPU (numpy) to avoid GPU memory pressure and pre-compiled
|
||
kernel mismatches (CUDA_ERROR_NO_BINARY_FOR_GPU on sm_89). The padded
|
||
array is transferred to GPU once, then each direction is processed and
|
||
results are streamed back to CPU.
|
||
|
||
Args:
|
||
dem: CPU numpy array — filled DEM (no NaN), shape (rows, cols).
|
||
rows, cols: dimensions.
|
||
res: resolution in m/px.
|
||
n_dirs: number of directions.
|
||
max_dist: max ray steps.
|
||
radii_m: list of radii in meters to record checkpoints.
|
||
If None, records only at max_dist.
|
||
|
||
Returns:
|
||
pos_angles: array of shape (n_dirs, n_radii, rows, cols) — max zenith angles
|
||
neg_angles: array of shape (n_dirs, n_radii, rows, cols) — max nadir angles
|
||
"""
|
||
angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
|
||
dx_dir = np.cos(angles)
|
||
dy_dir = np.sin(angles)
|
||
|
||
if radii_m is not None:
|
||
radii_steps = [min(int(r / res), max_dist) for r in radii_m]
|
||
n_radii = len(radii_m)
|
||
else:
|
||
radii_steps = [max_dist]
|
||
n_radii = 1
|
||
|
||
# Pad on CPU (numpy) — avoids GPU memory pressure and
|
||
# pre-compiled kernel issues (NO_BINARY_FOR_GPU on sm_89).
|
||
padded_np = np.pad(dem, max_dist, mode='constant', constant_values=np.nan)
|
||
|
||
# Transfer padded DEM to GPU for computation
|
||
padded = to_gpu(padded_np)
|
||
# GPU view of central region — reference elevation for ray-tracing
|
||
dem = padded[max_dist:max_dist+rows, max_dist:max_dist+cols]
|
||
# Free the CPU copy — we don't need it anymore
|
||
del padded_np
|
||
|
||
# Checkpoints triés par pas : (step, r_idx). Les snapshots sont pris quand
|
||
# le pas courant atteint le pas du checkpoint — les rayons ne dépassent
|
||
# donc pas le plus grand checkpoint demandé (équivalent au break d'avant).
|
||
checkpoints = sorted((radii_steps[r_idx], r_idx) for r_idx in range(n_radii))
|
||
last_step = checkpoints[-1][0]
|
||
|
||
# Process one direction at a time to limit GPU memory.
|
||
# Store results as flat CPU arrays — transfer back to GPU at the end.
|
||
pos_results = [None] * n_dirs
|
||
neg_results = [None] * n_dirs
|
||
|
||
for d_idx in range(n_dirs):
|
||
ddx, ddy = dx_dir[d_idx], dy_dir[d_idx]
|
||
|
||
# Pre-compute valid steps for this direction (jusqu'au dernier checkpoint)
|
||
valid_steps = []
|
||
for step in range(1, last_step + 1):
|
||
px = int(round(ddx * step))
|
||
py = int(round(ddy * step))
|
||
dist_m = math.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
|
||
if dist_m < res * 0.5:
|
||
continue
|
||
valid_steps.append((step, px, py, dist_m))
|
||
|
||
# Max cumulé des tangentes (float32 : moitié de VRAM vs float64)
|
||
running_pos = xp.zeros((rows, cols), dtype=np.float32)
|
||
running_neg = xp.zeros((rows, cols), dtype=np.float32)
|
||
snapshots = {}
|
||
cp_queue = list(checkpoints)
|
||
|
||
for step, px, py, dist_m in valid_steps:
|
||
# Slice from padded array, subtract original dem
|
||
view = padded[max_dist + py:max_dist + py + rows,
|
||
max_dist + px:max_dist + px + cols]
|
||
elev_diff = view - dem
|
||
del view # free slice reference
|
||
|
||
# Tangentes des angles (positive : terrain au-dessus, négative : en
|
||
# dessous). fmax propage le non-NaN : le bord de padding ne compte
|
||
# pas, comme avec l'ancien where(isnan) — en une seule opération.
|
||
running_pos = xp.fmax(running_pos,
|
||
xp.maximum(elev_diff, 0) / dist_m)
|
||
running_neg = xp.fmax(running_neg,
|
||
xp.maximum(-elev_diff, 0) / dist_m)
|
||
del elev_diff # free intermediate
|
||
|
||
# Snapshot du checkpoint atteint : conversion en angle UNE fois
|
||
while cp_queue and step >= cp_queue[0][0]:
|
||
_, r_idx = cp_queue.pop(0)
|
||
snapshots[r_idx] = (xp.arctan(running_pos),
|
||
xp.arctan(running_neg))
|
||
if not cp_queue:
|
||
break
|
||
|
||
# Checkpoints jamais atteints (steps invalides) : état final du balayage
|
||
while cp_queue:
|
||
_, r_idx = cp_queue.pop(0)
|
||
snapshots[r_idx] = (xp.arctan(running_pos),
|
||
xp.arctan(running_neg))
|
||
|
||
# Store results on CPU, free GPU memory before next direction
|
||
pos_results[d_idx] = to_cpu(xp.stack([snapshots[r][0] for r in range(n_radii)]))
|
||
neg_results[d_idx] = to_cpu(xp.stack([snapshots[r][1] for r in range(n_radii)]))
|
||
del running_pos, running_neg, snapshots
|
||
gpu_cleanup()
|
||
|
||
# Free the large padded array
|
||
del padded
|
||
gpu_cleanup()
|
||
|
||
# Reassemble into final arrays (on CPU to avoid GPU memory pressure)
|
||
pos_angles = np.array(pos_results)
|
||
neg_angles = np.array(neg_results)
|
||
|
||
return pos_angles, neg_angles
|
||
|
||
|
||
def _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m=None):
|
||
"""Ray-tracing avec repli CPU automatique si la VRAM est insuffisante.
|
||
|
||
Les dalles 0,2 m (5000×5000 px) multi-rayons peuvent dépasser la VRAM
|
||
disponible (GPU partagé avec d'autres services) : plutôt que d'abandonner
|
||
la visualisation, on désactive le GPU pour ce worker et on relance le
|
||
calcul sur CPU.
|
||
"""
|
||
try:
|
||
return _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
except Exception as e:
|
||
if _gpu_mod.is_gpu_active() and "out of memory" in str(e).lower():
|
||
logger.warning(" ⚠ VRAM insuffisante (ray-tracing) — repli CPU pour ce worker")
|
||
_gpu_mod.disable_gpu()
|
||
gpu_cleanup()
|
||
return _ray_trace_horizons_core(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
raise
|
||
|
||
|
||
# ============================================================
|
||
# Core terrain visualizations
|
||
# ============================================================
|
||
|
||
def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Generate multi-directional hillshade with contrast enhancement — GPU if available.
|
||
|
||
Combines 8-direction hillshade with slope shading for balanced illumination.
|
||
Applies percentile normalization and gamma correction to restore
|
||
contrast lost by averaging multiple azimuths.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Hillshade multidirectionnel{gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_hillshade_multi.tif"
|
||
used_gpu = _gpu_mod.HAS_GPU
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem = to_gpu(shared.dem_np)
|
||
dy = to_gpu(shared.dy) if _gpu_mod.HAS_GPU else shared.dy
|
||
dx = to_gpu(shared.dx) if _gpu_mod.HAS_GPU else shared.dx
|
||
slope = to_gpu(shared.slope_rad) if _gpu_mod.HAS_GPU else shared.slope_rad
|
||
aspect = xp.arctan2(dy, dx)
|
||
sin_slope = xp.sin(slope)
|
||
cos_slope = xp.cos(slope)
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
dem = to_gpu(dem_np)
|
||
dy, dx = xp.gradient(dem)
|
||
slope = xp.arctan(xp.sqrt(dx**2 + dy**2))
|
||
aspect = xp.arctan2(dy, dx)
|
||
sin_slope = xp.sin(slope)
|
||
cos_slope = xp.cos(slope)
|
||
|
||
# 8 azimuths for balanced illumination (eliminates directional bias)
|
||
azimuts = [0, 45, 90, 135, 180, 225, 270, 315]
|
||
altitude = 35 # Higher altitude for better micro-relief detection
|
||
hillshades = []
|
||
|
||
alt_rad = xp.radians(xp.array(altitude))
|
||
sin_alt = xp.sin(alt_rad)
|
||
cos_alt = xp.cos(alt_rad)
|
||
|
||
for az in azimuts:
|
||
az_rad = xp.radians(xp.array(az))
|
||
hs = sin_alt * sin_slope + cos_alt * cos_slope * xp.cos(az_rad - aspect)
|
||
hillshades.append(xp.clip(hs, 0, 1))
|
||
|
||
combined_hillshade = xp.mean(xp.array(hillshades), axis=0)
|
||
slope_shaded = cos_slope
|
||
combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
|
||
|
||
# Contrast enhancement: percentile stretch + gamma
|
||
combined_np = to_cpu(combined)
|
||
nan_mask = shared.nan_mask if shared else np.isnan(dem_np)
|
||
valid = combined_np[~nan_mask]
|
||
if len(valid) > 0:
|
||
p2, p98 = np.percentile(valid, 2), np.percentile(valid, 98)
|
||
if p98 - p2 > 0.01:
|
||
combined_np = np.clip((combined_np - p2) / (p98 - p2), 0, 1)
|
||
# Gamma correction to enhance shadows
|
||
gamma = 0.8
|
||
combined_np = np.power(combined_np, gamma)
|
||
|
||
_save_tif(output, combined_np.astype(np.float32), transform, crs, nan_mask=nan_mask)
|
||
logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur hillshade: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_slope(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Generate slope map (degrees) — GPU if available."""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Pente (Slope){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_slope.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
slope = shared.slope_deg
|
||
nan_mask = shared.nan_mask
|
||
if _gpu_mod.HAS_GPU:
|
||
slope = to_gpu(slope)
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
dem = to_gpu(dem_np)
|
||
dy, dx = xp.gradient(dem)
|
||
slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) * 180 / xp.pi
|
||
nan_mask = np.isnan(dem_np)
|
||
_save_tif(output, to_cpu(slope) if _gpu_mod.HAS_GPU else slope, transform, crs, nan_mask=nan_mask)
|
||
logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur slope: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Generate aspect (slope orientation) map — GPU if available.
|
||
|
||
0° = North, 90° = East, 180° = South, 270° = West.
|
||
Direction toward which the terrain descends.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Aspect (Orientation des pentes){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_aspect.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dy = shared.dy
|
||
dx = shared.dx
|
||
nan_mask = shared.nan_mask
|
||
if _gpu_mod.HAS_GPU:
|
||
dy = to_gpu(dy)
|
||
dx = to_gpu(dx)
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
dem = to_gpu(dem_np)
|
||
dy, dx = xp.gradient(dem)
|
||
nan_mask = np.isnan(dem_np)
|
||
aspect = xp.arctan2(dy, dx) * 180 / xp.pi
|
||
aspect = xp.mod(aspect, 360)
|
||
_save_tif(output, to_cpu(aspect) if _gpu_mod.HAS_GPU else aspect, transform, crs, nan_mask=nan_mask)
|
||
logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur aspect: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# GPU-accelerated visualizations
|
||
# ============================================================
|
||
|
||
def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Sky-View Factor - ray-tracing on 16 azimuths, multi-radius (GPU if available).
|
||
|
||
Traces rays in 16 directions at 3 radii (25, 50, 100m) and combines
|
||
with weights favoring medium range for archaeological feature detection.
|
||
SVF = (1/N) * sum(cos²(horizon_angle)). Valleys/crevices have low SVF,
|
||
ridges/peaks have high SVF.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Sky-View Factor (ray-tracing multi-rayon){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_svf.tif"
|
||
|
||
try:
|
||
dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
|
||
_prepare_dem_for_raycast(dem_file, shared, resolution)
|
||
|
||
radii_m = [25, 50, 100]
|
||
radius_weights = [0.3, 0.4, 0.3] # Medium range weighted more
|
||
max_dist = min(int(100 / res), 300)
|
||
n_dirs = 16
|
||
|
||
pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
|
||
# pos/neg are now numpy arrays (CPU) — combine on CPU
|
||
svf_combined = np.zeros((rows, cols), dtype=np.float32)
|
||
for r_idx in range(len(radii_m)):
|
||
horizon = np.maximum(pos_angles[:, r_idx], neg_angles[:, r_idx])
|
||
svf_r = np.mean(np.cos(horizon) ** 2, axis=0)
|
||
svf_combined += svf_r * radius_weights[r_idx]
|
||
|
||
svf_np = svf_combined
|
||
svf_np[nan_mask] = np.nan
|
||
_save_tif(output, svf_np, transform, crs)
|
||
logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur SVF: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None):
|
||
"""Positive/Negative Openness - multi-radius ray-tracing with std normalization.
|
||
|
||
Traces rays in 8 directions at 3 radii (25, 50, 100m).
|
||
Results are combined with equal weight across radii, then normalized
|
||
by standard deviation for cross-tile comparability.
|
||
"""
|
||
name = "positive_openness" if positive else "negative_openness"
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing multi-rayon){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_{name}.tif"
|
||
|
||
try:
|
||
dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
|
||
_prepare_dem_for_raycast(dem_file, shared, resolution)
|
||
|
||
radii_m = [25, 50, 100]
|
||
max_dist = min(int(100 / res), 300)
|
||
n_dirs = 8
|
||
|
||
pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
|
||
# Select positive or negative
|
||
if positive:
|
||
angles = pos_angles
|
||
else:
|
||
angles = neg_angles
|
||
|
||
# Mean across directions and radii (equal weight) — on CPU now
|
||
openness = np.mean(angles, axis=(0, 1))
|
||
openness_result = np.degrees(openness).astype(np.float32)
|
||
openness_result[nan_mask] = np.nan
|
||
|
||
# Z-score (écarts locaux en sigmas) : unités comparables entre tuiles,
|
||
# plage de rendu fixe → mosaïque de couleur homogène
|
||
valid = openness_result[~nan_mask]
|
||
if len(valid) > 0:
|
||
std_val = max(np.nanstd(valid), 0.01)
|
||
openness_result = (openness_result - np.nanmean(valid)) / std_val
|
||
|
||
_save_tif(output, openness_result, transform, crs)
|
||
logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur openness: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).
|
||
|
||
Scales adapt to resolution. Std normalization per scale.
|
||
Weighted combination favoring archaeologically relevant scales (5-25m).
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_mslrm.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
|
||
# Adaptive scales: finer at higher resolution
|
||
min_scale = max(2.0, resolution * 4)
|
||
# Archaeological scales: focus on 2-50m range.
|
||
# Small features (ditches, walls, post-holes) need 2-10m.
|
||
# Medium features (enclosures, roundhouses) need 10-25m.
|
||
# Large scales (50m+) are kept only for context with low weight.
|
||
candidate_scales = [2, 3, 5, 8, 10, 15, 25, 50]
|
||
sigmas = [s for s in candidate_scales if s >= min_scale]
|
||
|
||
# Weights: favor small-to-medium scales where archaeo features live
|
||
scale_weights = {
|
||
2: 1.5, 3: 1.8, 5: 2.0, 8: 1.8, 10: 1.5, 15: 1.3, 25: 1.0, 50: 0.5,
|
||
}
|
||
weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
|
||
|
||
logger.info(f" MSRM échelles: {sigmas}m")
|
||
lrm_stack = []
|
||
|
||
for sigma in sigmas:
|
||
sigma_px = sigma / resolution
|
||
if shared:
|
||
local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px)
|
||
else:
|
||
local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
|
||
lrm = dem_np - local_mean
|
||
lrm[nan_mask] = np.nan
|
||
# Std normalization
|
||
valid_lrm = lrm[~nan_mask]
|
||
lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
|
||
lrm = lrm / lrm_std
|
||
# Clip |z| to 3.0 to prevent large-scale outliers from drowning small features
|
||
lrm = np.clip(lrm, -3.0, 3.0)
|
||
lrm_stack.append(lrm.astype(np.float32))
|
||
|
||
# Weighted combination — preserve sign for RdBu_r colormap
|
||
# Positive = elevated (red), Negative = depression (blue)
|
||
lrm_array = np.array(lrm_stack)
|
||
weights_3d = weights[:, np.newaxis, np.newaxis]
|
||
with np.errstate(invalid='ignore', divide='ignore'):
|
||
with warnings.catch_warnings():
|
||
warnings.filterwarnings('ignore', message='Mean of empty slice')
|
||
# Signed RMS: magnitude from RMS, sign from weighted mean
|
||
signed_mean = np.nansum(lrm_array * weights_3d, axis=0) / np.sum(weights)
|
||
rms_magnitude = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
|
||
mslrm = np.sign(signed_mean) * rms_magnitude
|
||
mslrm[nan_mask] = np.nan
|
||
_save_tif(output, mslrm.astype(np.float32), transform, crs)
|
||
logger.info(f" ✓ MSRM terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur MSRM: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# SAILORE
|
||
# ============================================================
|
||
|
||
def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available).
|
||
|
||
Kernel size adapts to local slope: flat areas get larger kernels,
|
||
steep areas get smaller kernels. Scales adapt to resolution.
|
||
Reuses shared.lrm_15 when available to avoid recomputation.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → SAILORE (LRM adaptatif){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_sailore.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
slope_deg = shared.slope_deg
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
gy, gx = np.gradient(dem_np, resolution)
|
||
slope = np.arctan(np.sqrt(gx**2 + gy**2))
|
||
slope_deg = np.degrees(slope)
|
||
slope_deg[nan_mask] = np.nan
|
||
|
||
# Fixed physical scales (independent of resolution)
|
||
sigma_min_m = 2.0 # 2m — fine detail
|
||
sigma_max_m = 25.0 # 25m — broad relief
|
||
sigma_min = sigma_min_m / resolution
|
||
sigma_max = sigma_max_m / resolution
|
||
slope_norm = np.clip(slope_deg / 30.0, 0, 1)
|
||
|
||
# LRM fine (2m) — always compute
|
||
if shared:
|
||
lrm_fine = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_min)
|
||
else:
|
||
lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min)
|
||
lrm_fine[nan_mask] = np.nan
|
||
|
||
# LRM medium (13.5m) — reuse shared.lrm_15 (σ=15m) when available
|
||
sigma_mid = (sigma_min + sigma_max) / 2
|
||
if shared and abs(15.0 / resolution - sigma_mid) < 2.0 / resolution:
|
||
# shared.lrm_15 is close enough to medium scale
|
||
lrm_medium = shared.lrm_15.copy()
|
||
else:
|
||
if shared:
|
||
lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_mid)
|
||
else:
|
||
lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_mid)
|
||
lrm_medium[nan_mask] = np.nan
|
||
|
||
# LRM coarse (25m) — always compute
|
||
if shared:
|
||
lrm_coarse = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_max)
|
||
else:
|
||
lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max)
|
||
lrm_coarse[nan_mask] = np.nan
|
||
|
||
w_fine = slope_norm
|
||
w_medium = 1 - 2 * np.abs(slope_norm - 0.5)
|
||
w_coarse = 1 - slope_norm
|
||
w_total = w_fine + w_medium + w_coarse
|
||
w_total[w_total == 0] = 1
|
||
|
||
sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total
|
||
sailore[nan_mask] = np.nan
|
||
|
||
# Z-score (σ locales) : unités comparables entre tuiles, plage de
|
||
# rendu fixe ±3σ → mosaïque de couleur homogène
|
||
valid = sailore[~nan_mask]
|
||
if len(valid) > 0:
|
||
std_val = max(np.nanstd(valid), 0.01)
|
||
sailore = (sailore - np.nanmean(valid)) / std_val
|
||
|
||
_save_tif(output, sailore.astype(np.float32), transform, crs)
|
||
logger.info(f" ✓ SAILORE terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur SAILORE: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Roughness
|
||
# ============================================================
|
||
|
||
def _integral_sums(x):
|
||
"""Sommes intégrales 2D : S[i,j] = somme de x[0:i, 0:j] (float64).
|
||
|
||
Ligne/colonne 0 remplies de zéros — permet la somme d'une fenêtre
|
||
quelconque par 4 coins, y compris contre le bord (indices 0).
|
||
"""
|
||
rows, cols = x.shape
|
||
S = xp.zeros((rows + 1, cols + 1), dtype=np.float64)
|
||
S[1:, 1:] = xp.cumsum(xp.cumsum(x.astype(np.float64), axis=0), axis=1)
|
||
return S
|
||
|
||
|
||
def _box_std_from_integral(Sx, Sx2, size):
|
||
"""Écart-type local sur fenêtre size×size via sommes intégrales.
|
||
|
||
Coût indépendant de la taille de fenêtre (4 accès par pixel) — contre un
|
||
uniform_filter dont le coût croît avec la fenêtre (75 px à 0,2 m pour
|
||
l'échelle large). Aux bords, la fenêtre est tronquée et normalisée par le
|
||
nombre réel d'éléments (les dalles se recouvrent, le bord est sans effet
|
||
visuel).
|
||
"""
|
||
rows = Sx.shape[0] - 1
|
||
cols = Sx.shape[1] - 1
|
||
r = size // 2
|
||
iy0 = xp.maximum(xp.arange(rows) - r, 0)
|
||
iy1 = xp.minimum(xp.arange(rows) + r + 1, rows)
|
||
ix0 = xp.maximum(xp.arange(cols) - r, 0)
|
||
ix1 = xp.minimum(xp.arange(cols) + r + 1, cols)
|
||
# Nombre d'éléments réels de la fenêtre (tronquée aux bords)
|
||
counts = ((iy1 - iy0)[:, None] * (ix1 - ix0)[None, :]).astype(np.float64)
|
||
|
||
def box_sum(S):
|
||
return (S[iy1][:, ix1] - S[iy0][:, ix1]
|
||
- S[iy1][:, ix0] + S[iy0][:, ix0])
|
||
|
||
mean = box_sum(Sx) / counts
|
||
mean_sq = box_sum(Sx2) / counts
|
||
return xp.sqrt(xp.maximum(mean_sq - mean * mean, 0))
|
||
|
||
|
||
def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Surface roughness - multi-scale standard deviation (GPU-accelerated).
|
||
|
||
Combines fine (3m) and broad (15m) roughness for better detection
|
||
of archaeological features at multiple scales. Les écarts-types locaux
|
||
sont calculés par sommes intégrales : deux cumsum partagés entre les
|
||
deux échelles, extraction par 4 coins — coût constant quelle que soit
|
||
la fenêtre (un uniform_filter coûte proportionnellement à sa taille,
|
||
75 px à 0,2 m pour l'échelle large).
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Rugosité de surface{gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_roughness.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
if _gpu_mod.HAS_GPU:
|
||
filled = shared.filled_gpu
|
||
else:
|
||
filled = shared.filled
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
filled, _ = _fill_nans(dem_np)
|
||
if _gpu_mod.HAS_GPU:
|
||
filled = to_gpu(filled)
|
||
|
||
# Sommes intégrales partagées par les deux échelles (X et X²)
|
||
Sx = _integral_sums(filled)
|
||
Sx2 = _integral_sums(filled.astype(np.float64) ** 2)
|
||
|
||
fine_size = max(3, int(3 / resolution))
|
||
if fine_size % 2 == 0:
|
||
fine_size += 1
|
||
broad_size = max(3, int(15 / resolution))
|
||
if broad_size % 2 == 0:
|
||
broad_size += 1
|
||
|
||
roughness_fine = to_cpu(_box_std_from_integral(Sx, Sx2, fine_size))
|
||
roughness_broad = to_cpu(_box_std_from_integral(Sx, Sx2, broad_size))
|
||
del Sx, Sx2
|
||
gpu_cleanup()
|
||
roughness_fine[nan_mask] = np.nan
|
||
roughness_broad[nan_mask] = np.nan
|
||
|
||
# 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]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur rugosité: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Exposition des surfaces (Éclairage Solaire)
|
||
# ============================================================
|
||
|
||
def generate_solar(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Generate solar irradiance simulation.
|
||
|
||
Simulates morning sunlight (azimuth 90°, altitude 30°) to reveal
|
||
subtle topographic features through shadow effects.
|
||
"""
|
||
logger.info(" → Exposition des surfaces (Éclairage Solaire)...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_solar.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
nan_mask = shared.nan_mask
|
||
dx = shared.dx
|
||
dy = shared.dy
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
dem_filled, _ = _fill_nans(dem_np)
|
||
dy, dx = np.gradient(dem_filled, resolution, resolution)
|
||
|
||
# Solar parameters: morning sun (azimuth 90° = east, altitude 30°)
|
||
sun_azimuth = np.radians(90)
|
||
sun_altitude = np.radians(30)
|
||
|
||
# Aspect from gradient
|
||
aspect = np.degrees(np.arctan2(-dx, -dy))
|
||
aspect[aspect < 0] += 360
|
||
|
||
# Slope in radians
|
||
slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
|
||
|
||
# Solar irradiance calculation
|
||
irradiance = (np.sin(sun_altitude) * np.sin(slope_rad) +
|
||
np.cos(sun_altitude) * np.cos(slope_rad) *
|
||
np.cos(np.radians(aspect) - sun_azimuth))
|
||
|
||
# Clip to valid range [0, 1]
|
||
irradiance = np.clip(irradiance, 0, 1)
|
||
irradiance[nan_mask] = np.nan
|
||
|
||
_save_tif(output, irradiance.astype(np.float32), transform, crs)
|
||
logger.info(f" ✓ Exposition des surfaces terminée ({time.time()-t0:.1f}s)")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur exposition: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Wavelet (Mexican Hat)
|
||
# ============================================================
|
||
|
||
def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Mexican Hat wavelet multi-scale analysis (GPU if available).
|
||
|
||
Focused on small archaeological structures (paths, ditches, ramparts).
|
||
CWT 2D at scales [1, 2, 5, 10, 20, 50]m (0.5m added below 0.25m/px).
|
||
The 100m scale was dropped: it mostly responds to landforms (hills,
|
||
valleys), not to structures.
|
||
|
||
Large landforms are removed first by subtracting a Gaussian local-mean
|
||
trend (~35m). The residual is analyzed relative to its ~35m neighborhood,
|
||
so a ditch on a hilltop or slope does not stand out more than the same
|
||
ditch on flat ground (topographic-position independence). A Gaussian
|
||
smoothing preserves locally planar slopes, so slopes are removed too.
|
||
|
||
Uses robust per-scale normalization (MAD) and median-centered weighted RMS
|
||
combination with emphasis on small scales (1-10m).
|
||
The output is an index relative to the tile's own median level (≈ 1):
|
||
comparable from tile to tile, so a single fixed color stretch works.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Ondelette Mexican Hat multi-échelle{gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_wavelet.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
filled = shared.filled.astype(np.float64)
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
filled, _ = _fill_nans(dem_np.astype(np.float64))
|
||
|
||
min_scale = max(resolution * 2, 1.0)
|
||
# 100m retirée : elle répond surtout aux grandes formes du terrain,
|
||
# pas aux structures. Accent sur 1-10m (petites structures).
|
||
candidate_scales = [0.5, 1, 2, 5, 10, 20, 50]
|
||
scales = [s for s in candidate_scales if s >= min_scale]
|
||
|
||
scale_weights = {
|
||
0.5: 0.7, 1.0: 1.2, 2.0: 1.8, 5.0: 2.2,
|
||
10.0: 2.0, 20.0: 1.5, 50.0: 0.8,
|
||
}
|
||
weights = np.array([scale_weights.get(s, 1.0) for s in scales])
|
||
|
||
logger.info(f" Échelles CWT: {scales}m (résolution {resolution}m/px)")
|
||
|
||
from scipy.ndimage import gaussian_laplace, gaussian_filter
|
||
|
||
# Retrait des grands volumes (collines, vallées) : on soustrait une
|
||
# moyenne locale gaussienne avant la CWT. Un lissage gaussien préserve
|
||
# les pentes planes, donc le résidu est analysé par rapport à son
|
||
# voisinage ~35m : un fossé en sommet ou en flanc de colline ne
|
||
# ressort pas plus que le même fossé à plat.
|
||
# Fraction conservée pour une structure gaussienne de largeur σ_f :
|
||
# σ_t²/(σ_f²+σ_t²) → 10m : 92%, 20m : 75%, colline 150m : 5%.
|
||
# Mesuré sur MNT synthétique bruité : le contraste des petites
|
||
# structures est insensible à σ_t ; seul le fond sommet/plat varie
|
||
# (1.71 sans détendage → 1.21 à 35m).
|
||
detrend_sigma_m = 35.0
|
||
detrend_sigma_px = detrend_sigma_m / resolution
|
||
if _gpu_mod.HAS_GPU:
|
||
try:
|
||
from cupyx.scipy.ndimage import gaussian_filter as gpu_gaussian_filter
|
||
trend = to_cpu(gpu_gaussian_filter(to_gpu(filled), sigma=detrend_sigma_px))
|
||
except Exception:
|
||
trend = gaussian_filter(filled, sigma=detrend_sigma_px)
|
||
else:
|
||
trend = gaussian_filter(filled, sigma=detrend_sigma_px)
|
||
residual = filled - trend
|
||
del trend
|
||
logger.info(f" Retrait des grands volumes (tendance gaussienne {detrend_sigma_m:.0f}m)")
|
||
|
||
wavelet_stack = []
|
||
|
||
for scale_m in scales:
|
||
sigma_px = scale_m / resolution
|
||
if _gpu_mod.HAS_GPU:
|
||
try:
|
||
from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace
|
||
response = -gpu_gaussian_laplace(to_gpu(residual), sigma=sigma_px)
|
||
response = to_cpu(response)
|
||
except Exception:
|
||
response = -gaussian_laplace(residual, sigma=sigma_px)
|
||
else:
|
||
response = -gaussian_laplace(residual, sigma=sigma_px)
|
||
response[nan_mask] = np.nan
|
||
valid = response[~nan_mask]
|
||
# σ robuste (MAD) : le std classique est gonflé par les queues
|
||
# (structures marquées, bords de tuile) et varie fortement d'une
|
||
# tuile à l'autre — cause première des dominantes de couleur
|
||
# par tuile sur la carte.
|
||
mad = np.nanmedian(np.abs(valid - np.nanmedian(valid))) if len(valid) > 0 else 0.0
|
||
response = response / max(1.4826 * mad, 0.01)
|
||
wavelet_stack.append(response)
|
||
|
||
stack = np.array(wavelet_stack)
|
||
weights_3d = weights[:, np.newaxis, np.newaxis]
|
||
with np.errstate(invalid='ignore', divide='ignore'):
|
||
with warnings.catch_warnings():
|
||
warnings.filterwarnings('ignore', message='Mean of empty slice')
|
||
combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(weights))
|
||
combined[nan_mask] = np.nan
|
||
|
||
# Recentrage par la médiane de la tuile : la RMS devient un indice
|
||
# relatif au niveau moyen de la tuile (médiane = 1). La distribution
|
||
# est alors comparable d'une tuile à l'autre — condition pour un
|
||
# étirement couleur global fixe et homogène entre tuiles.
|
||
finite = combined[np.isfinite(combined)]
|
||
if finite.size:
|
||
combined = combined / max(float(np.median(finite)), 0.01)
|
||
combined[nan_mask] = np.nan
|
||
|
||
_save_tif(output, combined.astype(np.float32), transform, crs)
|
||
logger.info(f" ✓ Ondelette terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur ondelette: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Flow Accumulation helpers (module-level for numba caching)
|
||
# ============================================================
|
||
|
||
def _priority_flood(dem, nodata_mask):
|
||
"""Priority-flood algorithm for sink filling (Wang & Liu 2006).
|
||
|
||
O(n log n) compared to 50 iterations of minimum_filter.
|
||
Fills pits so water can flow downhill. NaN cells are treated as
|
||
closed (never pushed into the heap).
|
||
"""
|
||
result = _priority_flood_numba(dem, nodata_mask)
|
||
if result is not None:
|
||
return result
|
||
return _priority_flood_python(dem, nodata_mask)
|
||
|
||
|
||
def _priority_flood_python(dem, nodata_mask):
|
||
"""Pure-Python fallback for _priority_flood (used when numba is unavailable)."""
|
||
import heapq
|
||
|
||
rows, cols = dem.shape
|
||
filled = dem.copy()
|
||
closed = nodata_mask.copy()
|
||
open_queue = []
|
||
|
||
for r in range(rows):
|
||
for c in [0, cols - 1]:
|
||
if not closed[r, c]:
|
||
heapq.heappush(open_queue, (filled[r, c], r, c))
|
||
closed[r, c] = True
|
||
for c in range(1, cols - 1):
|
||
for r in [0, rows - 1]:
|
||
if not closed[r, c]:
|
||
heapq.heappush(open_queue, (filled[r, c], r, c))
|
||
closed[r, c] = True
|
||
|
||
dx8 = [1, 1, 0, -1, -1, -1, 0, 1]
|
||
dy8 = [0, 1, 1, 1, 0, -1, -1, -1]
|
||
|
||
while open_queue:
|
||
elev, r, c = heapq.heappop(open_queue)
|
||
for d in range(8):
|
||
nr, nc = r + dy8[d], c + dx8[d]
|
||
if 0 <= nr < rows and 0 <= nc < cols and not closed[nr, nc]:
|
||
if filled[nr, nc] < elev:
|
||
filled[nr, nc] = elev
|
||
closed[nr, nc] = True
|
||
heapq.heappush(open_queue, (filled[nr, nc], nr, nc))
|
||
|
||
return filled
|
||
|
||
|
||
def _priority_flood_numba(dem, nodata_mask):
|
||
"""JIT-compiled priority-flood via binary min-heap (~200x faster than Python).
|
||
|
||
Returns None if numba is unavailable (caller falls back to Python).
|
||
"""
|
||
try:
|
||
from numba import njit
|
||
except ImportError:
|
||
return None
|
||
|
||
@njit(cache=True)
|
||
def _flood(dem, nodata):
|
||
rows, cols = dem.shape
|
||
filled = dem.copy()
|
||
flat = filled.ravel()
|
||
closed = nodata.copy()
|
||
n = rows * cols
|
||
heap = np.empty(n, dtype=np.int64)
|
||
heap_size = 0
|
||
|
||
dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int8)
|
||
dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int8)
|
||
|
||
for r in range(rows):
|
||
for c in (0, cols - 1):
|
||
if not closed[r, c]:
|
||
heap[heap_size] = r * cols + c
|
||
heap_size += 1
|
||
closed[r, c] = True
|
||
for c in range(1, cols - 1):
|
||
for r in (0, rows - 1):
|
||
if not closed[r, c]:
|
||
heap[heap_size] = r * cols + c
|
||
heap_size += 1
|
||
closed[r, c] = True
|
||
|
||
while heap_size > 0:
|
||
cell = heap[0]
|
||
elev = flat[cell]
|
||
heap_size -= 1
|
||
if heap_size > 0:
|
||
heap[0] = heap[heap_size]
|
||
i = 0
|
||
while True:
|
||
l = 2 * i + 1
|
||
r = 2 * i + 2
|
||
smallest = i
|
||
if l < heap_size and flat[heap[l]] < flat[heap[smallest]]:
|
||
smallest = l
|
||
if r < heap_size and flat[heap[r]] < flat[heap[smallest]]:
|
||
smallest = r
|
||
if smallest == i:
|
||
break
|
||
heap[i], heap[smallest] = heap[smallest], heap[i]
|
||
i = smallest
|
||
|
||
r = cell // cols
|
||
c = cell % cols
|
||
for d in range(8):
|
||
nr = r + dy8[d]
|
||
nc = c + dx8[d]
|
||
if 0 <= nr < rows and 0 <= nc < cols and not closed[nr, nc]:
|
||
ncell = nr * cols + nc
|
||
if flat[ncell] < elev:
|
||
flat[ncell] = elev
|
||
closed[nr, nc] = True
|
||
heap[heap_size] = ncell
|
||
heap_size += 1
|
||
child = heap_size - 1
|
||
while child > 0:
|
||
parent = (child - 1) // 2
|
||
if flat[heap[child]] < flat[heap[parent]]:
|
||
heap[child], heap[parent] = heap[parent], heap[child]
|
||
child = parent
|
||
else:
|
||
break
|
||
|
||
return filled
|
||
|
||
return _flood(dem, nodata_mask)
|
||
|
||
|
||
def _d8_accumulate_numba(dem_filled, flow_dir, nodata_mask, rows, cols):
|
||
"""JIT-compiled D8 flow accumulation (top-down via elevation sort).
|
||
|
||
Uses numba for ~100x speedup over pure Python loop.
|
||
Falls back to pure Python if numba is unavailable.
|
||
"""
|
||
try:
|
||
from numba import njit
|
||
|
||
@njit(cache=True)
|
||
def _accumulate(dem, fdir, nodata, rows, cols):
|
||
dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int8)
|
||
dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int8)
|
||
|
||
flow_acc = np.ones((rows, cols), dtype=np.float32)
|
||
for r in range(rows):
|
||
for c in range(cols):
|
||
if nodata[r, c]:
|
||
flow_acc[r, c] = 0.0
|
||
|
||
# Sort cells by elevation descending (highest first)
|
||
n_cells = rows * cols
|
||
flat_dem = dem.ravel()
|
||
sort_idx = np.argsort(-flat_dem)
|
||
|
||
# Accumulate top-down (highest cell first)
|
||
for i in range(n_cells):
|
||
cell = sort_idx[i]
|
||
r = cell // cols
|
||
c = cell % cols
|
||
if nodata[r, c]:
|
||
continue
|
||
d = fdir[r, c]
|
||
if d < 0:
|
||
continue
|
||
nr = r + dy8[d]
|
||
nc = c + dx8[d]
|
||
if 0 <= nr < rows and 0 <= nc < cols and not nodata[nr, nc]:
|
||
flow_acc[nr, nc] += flow_acc[r, c]
|
||
|
||
return flow_acc
|
||
|
||
return _accumulate(dem_filled, flow_dir, nodata_mask, rows, cols)
|
||
|
||
except ImportError:
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Flow Accumulation
|
||
# ============================================================
|
||
|
||
def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Flow Accumulation — priority-flood sink filling + D8 accumulation.
|
||
|
||
Detects channels, ditches, and drainage paths by computing how many
|
||
upstream cells flow through each cell. Archaeological ditches and
|
||
natural drainage features both accumulate high flow values.
|
||
|
||
D8 direction is computed via vectorized numpy slicing.
|
||
Accumulation uses numba JIT (cached at module level) or pure Python fallback.
|
||
"""
|
||
logger.info(f" → Accumulation d'écoulement (flow accumulation)...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_flow_acc.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
filled = shared.filled
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
filled, _ = _fill_nans(dem_np)
|
||
|
||
rows, cols = dem_np.shape
|
||
|
||
# Sink filling — priority-flood (O(n log n), NaN-aware)
|
||
dem_filled = _priority_flood(filled, nan_mask)
|
||
|
||
logger.info(f" ✓ Sink filling terminé ({time.time()-t0:.1f}s)")
|
||
|
||
# D8 flow direction — vectorized via numpy slicing
|
||
dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int32)
|
||
dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int32)
|
||
dist8 = np.array([1.0, np.sqrt(2), 1.0, np.sqrt(2),
|
||
1.0, np.sqrt(2), 1.0, np.sqrt(2)])
|
||
|
||
flow_dir = np.full((rows, cols), -1, dtype=np.int8)
|
||
max_slope = np.zeros((rows, cols), dtype=np.float64)
|
||
|
||
padded = np.pad(dem_filled, 1, mode='constant',
|
||
constant_values=np.nanmax(dem_filled[~np.isnan(dem_filled)]) + 10000)
|
||
|
||
for d in range(8):
|
||
nx = 1 + dx8[d]
|
||
ny = 1 + dy8[d]
|
||
neighbor_elev = padded[ny:ny + rows, nx:nx + cols]
|
||
slope = (dem_filled - neighbor_elev) / (dist8[d] * resolution)
|
||
slope[nan_mask] = -1
|
||
better = slope > max_slope
|
||
flow_dir[better] = d
|
||
max_slope[better] = slope[better]
|
||
|
||
logger.info(f" ✓ Direction D8 terminée ({time.time()-t0:.1f}s)")
|
||
|
||
# D8 accumulation — numba JIT (module-level cache) or pure Python fallback
|
||
result = _d8_accumulate_numba(dem_filled, flow_dir, nan_mask.astype(np.bool_), rows, cols)
|
||
|
||
if result is not None:
|
||
flow_acc = result
|
||
logger.info(f" Accumulation D8 via numba")
|
||
else:
|
||
# Pure Python fallback
|
||
logger.info(f" Accumulation D8 via Python (installez numba pour accélérer)")
|
||
flat_dem = dem_filled[~nan_mask].flatten()
|
||
valid_indices = np.where(~nan_mask.flatten())[0]
|
||
sort_order = valid_indices[np.argsort(-flat_dem)]
|
||
|
||
flow_acc = np.ones((rows, cols), dtype=np.float32)
|
||
flow_acc[nan_mask] = 0
|
||
|
||
for idx in sort_order:
|
||
r, c = divmod(idx, cols)
|
||
d = flow_dir[r, c]
|
||
if d < 0:
|
||
continue
|
||
nr, nc = r + dy8[d], c + dx8[d]
|
||
if 0 <= nr < rows and 0 <= nc < cols and not nan_mask[nr, nc]:
|
||
flow_acc[nr, nc] += flow_acc[r, c]
|
||
|
||
logger.info(f" ✓ D8 accumulation terminé ({time.time()-t0:.1f}s)")
|
||
|
||
# Log transform
|
||
flow_result = np.log1p(flow_acc)
|
||
flow_result[nan_mask] = np.nan
|
||
|
||
_save_tif(output, flow_result, transform, crs)
|
||
logger.info(f" ✓ Flow accumulation terminé ({time.time()-t0:.1f}s)")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur flow accumulation: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Anomaly Mask — automatic threshold detection
|
||
# ============================================================
|
||
|
||
def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, n_sigma=2.0):
|
||
"""Composite anomaly mask — automatic threshold detection (GPU if available).
|
||
|
||
Reads pre-computed visualization layers (MSRM, SVF, Wavelet, Openness Neg,
|
||
Roughness), normalizes each to z-scores, and combines them into a composite
|
||
anomaly score. Pixels beyond `n_sigma` standard deviations of the local mean
|
||
are flagged as suspicious.
|
||
|
||
The output is a continuous score (0–1) where:
|
||
- 0 = no anomaly (flat/natural terrain)
|
||
- 1 = high anomaly (potential archaeological structure)
|
||
|
||
This mask is directly usable in GIS for polygon extraction and field survey
|
||
planning.
|
||
|
||
Args:
|
||
n_sigma: Number of standard deviations for the anomaly threshold.
|
||
Lower = more sensitive (more false positives).
|
||
Default 2.0 (good balance for archaeological detection).
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Détection automatique d'anomalies (seuil {n_sigma}σ){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_anomaly.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
nan_mask = shared.nan_mask
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
|
||
rows, cols = nan_mask.shape
|
||
|
||
# Collect available visualization layers from disk
|
||
# Each is loaded, converted to |z-score|, and contributes to a weighted sum.
|
||
# The weighted sum of |z| acts as a "vote": pixels where multiple layers
|
||
# show anomalies get higher scores than pixels where only one layer fires.
|
||
layer_configs = [
|
||
# (filename_pattern, weight)
|
||
("mslrm", 2.5), # Multi-scale relief — strongest signal
|
||
("negative_openness", 2.0), # Fossés, dolines
|
||
("roughness", 1.8), # Surface irregularity
|
||
("wavelet", 1.5), # Circular + linear structures
|
||
("svf", 1.3), # Sky-view depressions
|
||
("flow_acc", 1.2), # Drainage channels / ditches
|
||
("positive_openness", 0.8), # Surélevations
|
||
]
|
||
|
||
layers = []
|
||
for pattern, weight in layer_configs:
|
||
layer_path = vis_dir / f"{basename}_{pattern}.tif"
|
||
if not layer_path.exists():
|
||
layer_path = vis_dir / f"{basename}_negative_openness.tif" if "neg" in pattern else None
|
||
if layer_path is None or not layer_path.exists():
|
||
continue
|
||
try:
|
||
with rasterio.open(layer_path) as src:
|
||
data = src.read(1).astype(np.float64)
|
||
# Absolute z-score (captures both positive and negative deviations)
|
||
valid = data[~nan_mask]
|
||
if len(valid) == 0:
|
||
continue
|
||
mean_val = np.nanmean(valid)
|
||
std_val = max(np.nanstd(valid), 0.01)
|
||
zscore = np.abs(data - mean_val) / std_val
|
||
zscore[nan_mask] = 0.0
|
||
layers.append((zscore, weight))
|
||
except Exception as e:
|
||
logger.debug(f" Couche {pattern} non disponible: {e}")
|
||
continue
|
||
|
||
if not layers:
|
||
# Fallback: use MSRM + roughness computed on the fly
|
||
logger.info(" Aucune couche trouvée — calcul MSRM + rugosité en direct...")
|
||
dem_np_safe = shared.dem_np if shared else dem_np
|
||
|
||
# Quick MSRM (single scale 10m for speed)
|
||
sigma_px = max(5, 10.0 / resolution)
|
||
local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px) if shared else \
|
||
_filter_nanaware(dem_np_safe, xp_gaussian_filter, sigma=sigma_px)
|
||
quick_mslrm = np.abs(dem_np_safe - local_mean)
|
||
quick_mslrm[nan_mask] = np.nan
|
||
valid = quick_mslrm[~nan_mask]
|
||
std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
|
||
quick_mslrm = quick_mslrm / std_val
|
||
layers.append((quick_mslrm, 2.5))
|
||
|
||
# Quick roughness
|
||
fine_size = max(3, int(3 / resolution))
|
||
if fine_size % 2 == 0:
|
||
fine_size += 1
|
||
if shared:
|
||
fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
|
||
fine_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
|
||
else:
|
||
fine_mean = _filter_nanaware(dem_np_safe.astype(np.float64), xp_uniform_filter, size=fine_size)
|
||
fine_mean_sq = _filter_nanaware(dem_np_safe.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
|
||
roughness = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0))
|
||
roughness[nan_mask] = np.nan
|
||
valid_r = roughness[~nan_mask]
|
||
std_val_r = max(np.nanstd(valid_r), 0.01) if len(valid_r) > 0 else 0.01
|
||
roughness = roughness / std_val_r
|
||
layers.append((roughness, 1.5))
|
||
|
||
# Weighted SUM of |z-scores| (not RMS — each layer votes independently)
|
||
combined = np.zeros((rows, cols), dtype=np.float64)
|
||
total_weight = 0.0
|
||
for layer_data, weight in layers:
|
||
combined += layer_data * weight
|
||
total_weight += weight
|
||
|
||
if total_weight > 0:
|
||
combined = combined / total_weight
|
||
|
||
# Adaptive threshold: suppress pixels below the (100 - n_sigma*10)th percentile.
|
||
# With n_sigma=2.0 → 80th percentile: keep only the top 20% of signal.
|
||
# This adapts to each tile's terrain instead of a fixed z-score cutoff.
|
||
threshold_pct = max(50, min(95, 100 - n_sigma * 10))
|
||
threshold_val = np.percentile(combined[~nan_mask], threshold_pct)
|
||
combined = np.clip(combined - threshold_val, 0, None)
|
||
|
||
# Rescale survivors to 0–1
|
||
above_thresh = combined[combined > 0]
|
||
if len(above_thresh) > 0:
|
||
p95 = np.percentile(above_thresh, 95)
|
||
if p95 > 0:
|
||
combined = np.clip(combined / p95, 0, 1)
|
||
|
||
combined[nan_mask] = np.nan
|
||
|
||
_save_tif(output, combined.astype(np.float32), transform, crs)
|
||
n_anomaly = int(np.sum(combined > 0.1)) if np.any(combined > 0) else 0
|
||
pct_anomaly = n_anomaly / max(np.sum(~nan_mask), 1) * 100
|
||
logger.info(f" ✓ Détection anomalies terminée ({time.time()-t0:.1f}s) — "
|
||
f"{pct_anomaly:.1f}% de la zone ({n_anomaly} px) au-delà de {n_sigma}σ")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur détection anomalies: {e}", exc_info=True)
|
||
return None
|