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