0d3c3ba8bc
Nouvelle visualisation relief_oriente : une image RGB unique qui fusionne l'openness positive locale (MNT détendancé σ 10 m, rayons 5/10/20 m, 16 directions) portée par la clarté CIELAB et l'orientation des pentes portée par la teinte. Échelle log fixe et support de 40 m sous la bande de raccord de 100 m : dalles jointives. Calcul sur grille décimée à 0,8 m, noyau dédié (CuPy RawKernel, numba parallèle, repli numpy) et colorisation par table L* × teinte : ~8 s par dalle sur CPU au lieu de ~50 s. L'openness positive et négative est normalisée par des références figées mesurées sur 15 dalles au lieu d'un z-score par dalle, qui rendait l'échelle de couleur non jointive. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2006 lines
82 KiB
Python
2006 lines
82 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
|
||
# 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
|