1271 lines
50 KiB
Python
1271 lines
50 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 HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, xp_minimum_filter, xp_maximum_filter, gpu_cleanup
|
||
from . import gpu as _gpu_mod
|
||
|
||
logger = logging.getLogger("lidar")
|
||
|
||
# CuPy module reference — lazily imported on first GPU use.
|
||
# If disable_gpu() is called at runtime, HAS_GPU becomes False
|
||
# and xp delegates to numpy instead.
|
||
_cp = None
|
||
|
||
|
||
class _XPProxy:
|
||
"""Proxy that delegates array operations to cupy or numpy.
|
||
|
||
Checks HAS_GPU on every attribute access so that disable_gpu()
|
||
(called on CUDA errors) takes effect immediately, without needing
|
||
to change every call site in visualizations.py.
|
||
"""
|
||
|
||
def __getattr__(self, name):
|
||
global _cp
|
||
from . import gpu as _gpu_mod
|
||
if _gpu_mod.HAS_GPU:
|
||
if _cp is None:
|
||
try:
|
||
import cupy
|
||
_cp = cupy
|
||
except ImportError:
|
||
pass
|
||
if _cp is not None:
|
||
return getattr(_cp, name)
|
||
return getattr(np, name)
|
||
|
||
|
||
xp = _XPProxy()
|
||
|
||
|
||
class SharedDEM:
|
||
"""Pre-computed DEM data shared across all visualizations.
|
||
|
||
Reads the DEM once and lazily computes on first access:
|
||
- NaN mask and filled DEM (avoids 20+ calls to _fill_nans)
|
||
- Gradient components (shared by hillshade, slope)
|
||
- 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.
|
||
"""
|
||
from scipy.interpolate import NearestNDInterpolator
|
||
nan_mask = np.isnan(arr)
|
||
if not np.any(nan_mask):
|
||
return arr, nan_mask
|
||
valid = ~nan_mask
|
||
y_coords, x_coords = np.where(valid)
|
||
if len(y_coords) == 0:
|
||
return np.zeros_like(arr), nan_mask
|
||
z_values = arr[valid]
|
||
interp = NearestNDInterpolator(
|
||
np.column_stack((y_coords, x_coords)), z_values
|
||
)
|
||
y_missing, x_missing = np.where(nan_mask)
|
||
filled = arr.copy()
|
||
filled[y_missing, x_missing] = interp(y_missing, x_missing)
|
||
return filled, nan_mask
|
||
|
||
|
||
def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
|
||
"""Apply a filter to an array while preserving NaN zones.
|
||
|
||
1. Fill NaN with nearest-neighbor interpolation
|
||
2. Apply the filter
|
||
3. Restore original NaN mask on the result
|
||
|
||
Args:
|
||
arr: Input array (numpy or cupy).
|
||
filter_func: Function that takes (array, *args, **kwargs) and returns filtered array.
|
||
use_gpu: If True, apply filter on GPU (send filled array to GPU first).
|
||
Returns:
|
||
Filtered array with original NaN positions preserved.
|
||
"""
|
||
is_gpu_arr = _gpu_mod.HAS_GPU and _cp is not None and isinstance(arr, _cp.ndarray)
|
||
arr_np = to_cpu(arr) if is_gpu_arr else arr
|
||
|
||
filled, nan_mask = _fill_nans(arr_np)
|
||
|
||
if use_gpu and _gpu_mod.HAS_GPU:
|
||
filled_gpu = to_gpu(filled)
|
||
result_gpu = filter_func(filled_gpu, *args, **kwargs)
|
||
result = to_cpu(result_gpu)
|
||
gpu_cleanup()
|
||
else:
|
||
result = filter_func(filled, *args, **kwargs)
|
||
|
||
result[nan_mask] = np.nan
|
||
return result
|
||
|
||
|
||
# ============================================================
|
||
# Shared ray-tracing core
|
||
# ============================================================
|
||
|
||
def _prepare_dem_for_raycast(dem_file, shared, resolution):
|
||
"""Load DEM and prepare padded array for ray-tracing.
|
||
|
||
Returns (dem_gpu_or_cpu, dem_np, rows, cols, res, nan_mask,
|
||
transform, crs, padded) ready for ray-tracing.
|
||
"""
|
||
if shared:
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
filled, _ = _fill_nans(dem_np)
|
||
dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled
|
||
res = resolution
|
||
rows, cols = dem_np.shape
|
||
return dem, dem_np, rows, cols, res, nan_mask, transform, crs
|
||
|
||
|
||
def _ray_trace_horizons(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.
|
||
|
||
Args:
|
||
dem: GPU or CPU filled DEM array (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
|
||
|
||
padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
|
||
|
||
# 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
|
||
valid_steps = []
|
||
for step in range(1, max_dist + 1):
|
||
px = int(round(ddx * step))
|
||
py = int(round(ddy * step))
|
||
dist_m = math.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
|
||
if dist_m < res * 0.5:
|
||
continue
|
||
valid_steps.append((step, px, py, dist_m))
|
||
|
||
# Running max angles per radius
|
||
running_pos = xp.zeros((n_radii, rows, cols))
|
||
running_neg = xp.zeros((n_radii, rows, cols))
|
||
# Track which radius checkpoints have been passed
|
||
radii_remaining = set(range(n_radii))
|
||
|
||
for step, px, py, dist_m in valid_steps:
|
||
elev_diff = padded[max_dist + py:max_dist + py + rows,
|
||
max_dist + px:max_dist + px + cols] - dem
|
||
|
||
# Positive: angle to terrain above viewer
|
||
pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m)
|
||
# Negative: angle to terrain below viewer
|
||
neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m)
|
||
|
||
# Update running max for all radius checkpoints still active
|
||
for r_idx in radii_remaining:
|
||
pos_angle_safe = xp.nan_to_num(pos_angle, nan=0)
|
||
neg_angle_safe = xp.nan_to_num(neg_angle, nan=0)
|
||
running_pos[r_idx] = xp.where(xp.isnan(pos_angle), running_pos[r_idx],
|
||
xp.maximum(running_pos[r_idx], pos_angle_safe))
|
||
running_neg[r_idx] = xp.where(xp.isnan(neg_angle), running_neg[r_idx],
|
||
xp.maximum(running_neg[r_idx], neg_angle_safe))
|
||
|
||
# Check which radii have been passed
|
||
new_remaining = set()
|
||
for r_idx in radii_remaining:
|
||
if step < radii_steps[r_idx]:
|
||
new_remaining.add(r_idx)
|
||
radii_remaining = new_remaining
|
||
if not radii_remaining:
|
||
break
|
||
|
||
# Store results on CPU, free GPU memory before next direction
|
||
pos_results[d_idx] = to_cpu(running_pos)
|
||
neg_results[d_idx] = to_cpu(running_neg)
|
||
del running_pos, running_neg
|
||
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
|
||
|
||
|
||
# ============================================================
|
||
# Core terrain visualizations
|
||
# ============================================================
|
||
|
||
def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Generate multi-directional hillshade with contrast enhancement — GPU if available.
|
||
|
||
Combines 8-direction hillshade with slope shading for balanced illumination.
|
||
Applies percentile normalization and gamma correction to restore
|
||
contrast lost by averaging multiple azimuths.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Hillshade multidirectionnel{gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_hillshade_multi.tif"
|
||
used_gpu = _gpu_mod.HAS_GPU
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem = to_gpu(shared.dem_np)
|
||
dy = to_gpu(shared.dy) if _gpu_mod.HAS_GPU else shared.dy
|
||
dx = to_gpu(shared.dx) if _gpu_mod.HAS_GPU else shared.dx
|
||
slope = to_gpu(shared.slope_rad) if _gpu_mod.HAS_GPU else shared.slope_rad
|
||
aspect = xp.arctan2(dy, dx)
|
||
sin_slope = xp.sin(slope)
|
||
cos_slope = xp.cos(slope)
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
dem = to_gpu(dem_np)
|
||
dy, dx = xp.gradient(dem)
|
||
slope = xp.arctan(xp.sqrt(dx**2 + dy**2))
|
||
aspect = xp.arctan2(dy, dx)
|
||
sin_slope = xp.sin(slope)
|
||
cos_slope = xp.cos(slope)
|
||
|
||
# 8 azimuths for balanced illumination (eliminates directional bias)
|
||
azimuts = [0, 45, 90, 135, 180, 225, 270, 315]
|
||
altitude = 35 # Higher altitude for better micro-relief detection
|
||
hillshades = []
|
||
|
||
alt_rad = xp.radians(xp.array(altitude))
|
||
sin_alt = xp.sin(alt_rad)
|
||
cos_alt = xp.cos(alt_rad)
|
||
|
||
for az in azimuts:
|
||
az_rad = xp.radians(xp.array(az))
|
||
hs = sin_alt * sin_slope + cos_alt * cos_slope * xp.cos(az_rad - aspect)
|
||
hillshades.append(xp.clip(hs, 0, 1))
|
||
|
||
combined_hillshade = xp.mean(xp.array(hillshades), axis=0)
|
||
slope_shaded = cos_slope
|
||
combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
|
||
|
||
# Contrast enhancement: percentile stretch + gamma
|
||
combined_np = to_cpu(combined)
|
||
nan_mask = shared.nan_mask if shared else np.isnan(dem_np)
|
||
valid = combined_np[~nan_mask]
|
||
if len(valid) > 0:
|
||
p2, p98 = np.percentile(valid, 2), np.percentile(valid, 98)
|
||
if p98 - p2 > 0.01:
|
||
combined_np = np.clip((combined_np - p2) / (p98 - p2), 0, 1)
|
||
# Gamma correction to enhance shadows
|
||
gamma = 0.8
|
||
combined_np = np.power(combined_np, gamma)
|
||
|
||
_save_tif(output, combined_np.astype(np.float32), transform, crs, nan_mask=nan_mask)
|
||
logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur hillshade: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_slope(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Generate slope map (degrees) — GPU if available."""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Pente (Slope){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_slope.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
slope = shared.slope_deg
|
||
nan_mask = shared.nan_mask
|
||
if _gpu_mod.HAS_GPU:
|
||
slope = to_gpu(slope)
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
dem = to_gpu(dem_np)
|
||
dy, dx = xp.gradient(dem)
|
||
slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) * 180 / xp.pi
|
||
nan_mask = np.isnan(dem_np)
|
||
_save_tif(output, to_cpu(slope) if _gpu_mod.HAS_GPU else slope, transform, crs, nan_mask=nan_mask)
|
||
logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur slope: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# GPU-accelerated visualizations
|
||
# ============================================================
|
||
|
||
def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Sky-View Factor - ray-tracing on 16 azimuths, multi-radius (GPU if available).
|
||
|
||
Traces rays in 16 directions at 3 radii (25, 50, 100m) and combines
|
||
with weights favoring medium range for archaeological feature detection.
|
||
SVF = (1/N) * sum(cos²(horizon_angle)). Valleys/crevices have low SVF,
|
||
ridges/peaks have high SVF.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Sky-View Factor (ray-tracing multi-rayon){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_svf.tif"
|
||
|
||
try:
|
||
dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
|
||
_prepare_dem_for_raycast(dem_file, shared, resolution)
|
||
|
||
radii_m = [25, 50, 100]
|
||
radius_weights = [0.3, 0.4, 0.3] # Medium range weighted more
|
||
max_dist = min(int(100 / res), 300)
|
||
n_dirs = 16
|
||
|
||
pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
|
||
# pos/neg are now numpy arrays (CPU) — combine on CPU
|
||
svf_combined = np.zeros((rows, cols), dtype=np.float32)
|
||
for r_idx in range(len(radii_m)):
|
||
horizon = np.maximum(pos_angles[:, r_idx], neg_angles[:, r_idx])
|
||
svf_r = np.mean(np.cos(horizon) ** 2, axis=0)
|
||
svf_combined += svf_r * radius_weights[r_idx]
|
||
|
||
svf_np = svf_combined
|
||
svf_np[nan_mask] = np.nan
|
||
_save_tif(output, svf_np, transform, crs)
|
||
logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur SVF: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None):
|
||
"""Positive/Negative Openness - multi-radius ray-tracing with std normalization.
|
||
|
||
Traces rays in 8 directions at 3 radii (25, 50, 100m).
|
||
Results are combined with equal weight across radii, then normalized
|
||
by standard deviation for cross-tile comparability.
|
||
"""
|
||
name = "positive_openness" if positive else "negative_openness"
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing multi-rayon){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_{name}.tif"
|
||
|
||
try:
|
||
dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
|
||
_prepare_dem_for_raycast(dem_file, shared, resolution)
|
||
|
||
radii_m = [25, 50, 100]
|
||
max_dist = min(int(100 / res), 300)
|
||
n_dirs = 8
|
||
|
||
pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
|
||
# Select positive or negative
|
||
if positive:
|
||
angles = pos_angles
|
||
else:
|
||
angles = neg_angles
|
||
|
||
# Mean across directions and radii (equal weight) — on CPU now
|
||
openness = np.mean(angles, axis=(0, 1))
|
||
openness_result = np.degrees(openness).astype(np.float32)
|
||
openness_result[nan_mask] = np.nan
|
||
|
||
# Std normalization for cross-tile comparability
|
||
valid = openness_result[~nan_mask]
|
||
if len(valid) > 0:
|
||
std_val = max(np.nanstd(valid), 0.01)
|
||
openness_result = openness_result / std_val
|
||
|
||
_save_tif(output, openness_result, transform, crs)
|
||
logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur openness: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).
|
||
|
||
Scales adapt to resolution. Std normalization per scale.
|
||
Weighted combination favoring archaeologically relevant scales (5-25m).
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_mslrm.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
|
||
# Adaptive scales: finer at higher resolution
|
||
min_scale = max(2.0, resolution * 4)
|
||
candidate_scales = [2, 5, 10, 20, 50, 100, 200]
|
||
sigmas = [s for s in candidate_scales if s >= min_scale]
|
||
|
||
# Archaeological weights: favor 5-25m range (ditches, enclosures, tumulus)
|
||
scale_weights = {
|
||
2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6, 200: 0.4,
|
||
}
|
||
weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
|
||
|
||
logger.info(f" MSRM échelles: {sigmas}m")
|
||
lrm_stack = []
|
||
|
||
for sigma in sigmas:
|
||
sigma_px = sigma / resolution
|
||
if shared:
|
||
local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px)
|
||
else:
|
||
local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
|
||
lrm = dem_np - local_mean
|
||
lrm[nan_mask] = np.nan
|
||
# Std normalization: x / std — preserves sign and contrast better than z-score
|
||
valid_lrm = lrm[~nan_mask]
|
||
lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
|
||
lrm = lrm / lrm_std
|
||
lrm_stack.append(lrm.astype(np.float32))
|
||
|
||
# Weighted combination — preserve sign for RdBu_r colormap
|
||
# Positive = elevated (red), Negative = depression (blue)
|
||
lrm_array = np.array(lrm_stack)
|
||
weights_3d = weights[:, np.newaxis, np.newaxis]
|
||
with np.errstate(invalid='ignore', divide='ignore'):
|
||
with warnings.catch_warnings():
|
||
warnings.filterwarnings('ignore', message='Mean of empty slice')
|
||
# Signed RMS: magnitude from RMS, sign from weighted mean
|
||
signed_mean = np.nansum(lrm_array * weights_3d, axis=0) / np.sum(weights)
|
||
rms_magnitude = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
|
||
mslrm = np.sign(signed_mean) * rms_magnitude
|
||
mslrm[nan_mask] = np.nan
|
||
_save_tif(output, mslrm.astype(np.float32), transform, crs)
|
||
logger.info(f" ✓ MSRM terminé ({time.time()-t0:.1f}s){' [GPU]' 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
|
||
|
||
_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 generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Surface roughness - multi-scale standard deviation (GPU-accelerated).
|
||
|
||
Combines fine (3m) and broad (15m) roughness for better detection
|
||
of archaeological features at multiple scales.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Rugosité de surface{gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_roughness.tif"
|
||
|
||
try:
|
||
if shared:
|
||
transform = shared.transform
|
||
crs = shared.crs
|
||
dem_np = shared.dem_np
|
||
nan_mask = shared.nan_mask
|
||
else:
|
||
dem_np, transform, crs = _read_dem(dem_file)
|
||
nan_mask = np.isnan(dem_np)
|
||
|
||
# Fine roughness (3m window)
|
||
fine_size = max(3, int(3 / resolution))
|
||
if fine_size % 2 == 0:
|
||
fine_size += 1
|
||
|
||
if shared:
|
||
fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
|
||
fine_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
|
||
fine_mean_sq[shared.nan_mask] = np.nan
|
||
else:
|
||
fine_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=fine_size)
|
||
fine_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
|
||
roughness_fine = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0))
|
||
roughness_fine[nan_mask] = np.nan
|
||
|
||
# Broad roughness (15m window)
|
||
broad_size = max(3, int(15 / resolution))
|
||
if broad_size % 2 == 0:
|
||
broad_size += 1
|
||
|
||
if shared:
|
||
broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size)
|
||
broad_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
|
||
broad_mean_sq[shared.nan_mask] = np.nan
|
||
else:
|
||
broad_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=broad_size)
|
||
broad_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
|
||
roughness_broad = np.sqrt(np.maximum(broad_mean_sq - broad_mean * broad_mean, 0))
|
||
roughness_broad[nan_mask] = np.nan
|
||
|
||
# Std normalization per scale then weighted combination
|
||
fine_valid = roughness_fine[~nan_mask]
|
||
broad_valid = roughness_broad[~nan_mask]
|
||
fine_std = max(np.nanstd(fine_valid), 0.01) if len(fine_valid) > 0 else 0.01
|
||
broad_std = max(np.nanstd(broad_valid), 0.01) if len(broad_valid) > 0 else 0.01
|
||
|
||
roughness = 0.7 * roughness_fine / fine_std + 0.3 * roughness_broad / broad_std
|
||
roughness[nan_mask] = np.nan
|
||
|
||
_save_tif(output, roughness, transform, crs)
|
||
logger.info(f" ✓ Rugosité terminée ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur rugosité: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Wavelet (Mexican Hat + Directional Gabor)
|
||
# ============================================================
|
||
|
||
def _gabor_kernel_2d(size, sigma, wavelength, theta):
|
||
"""Create a 2D Gabor kernel.
|
||
|
||
Args:
|
||
size: kernel size (odd integer)
|
||
sigma: standard deviation
|
||
wavelength: wavelength of sinusoid
|
||
theta: orientation angle in radians (0 = horizontal)
|
||
"""
|
||
center = size // 2
|
||
y, x = np.ogrid[-center:center+1, -center:center+1]
|
||
# Rotate coordinates
|
||
x_theta = x * np.cos(theta) + y * np.sin(theta)
|
||
y_theta = -x * np.sin(theta) + y * np.cos(theta)
|
||
sigma_sq = 2 * sigma * sigma
|
||
kernel = np.exp(-(x_theta**2 + y_theta**2) / sigma_sq)
|
||
kernel *= np.cos(2 * np.pi * x_theta / wavelength)
|
||
return kernel
|
||
|
||
|
||
def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Multi-scale wavelet analysis: Mexican Hat + Directional Gabor (GPU if available).
|
||
|
||
Mexican Hat (radial): detects circular features (tumulus, enclos ronds).
|
||
Gabor (directional): detects linear features (chemins, murs, fossés).
|
||
4 Gabor orientations (0°, 45°, 90°, 135°) at key archaeological scales.
|
||
Both combined with RMS for orientation-invariant detection.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Ondelette Mexican Hat + Gabor directionnelle{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))
|
||
|
||
# --- Mexican Hat scales ---
|
||
min_scale = max(resolution * 2, 1.0)
|
||
candidate_scales = [0.5, 1, 2, 5, 10, 20, 50, 100]
|
||
mex_scales = [s for s in candidate_scales if s >= min_scale]
|
||
|
||
mex_weights_map = {
|
||
0.5: 0.6, 1.0: 0.8, 2.0: 1.5, 5.0: 2.0,
|
||
10.0: 1.8, 20.0: 1.5, 50.0: 1.0, 100.0: 0.6,
|
||
}
|
||
mex_weights = np.array([mex_weights_map.get(s, 1.0) for s in mex_scales])
|
||
|
||
logger.info(f" Échelles CWT: {mex_scales}m (résolution {resolution}m/px)")
|
||
|
||
from scipy.ndimage import gaussian_laplace, convolve
|
||
|
||
wavelet_stack = []
|
||
|
||
# Mexican Hat (radial) — multi-scale
|
||
for scale_m in mex_scales:
|
||
sigma_px = scale_m / resolution
|
||
if _gpu_mod.HAS_GPU:
|
||
try:
|
||
from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace
|
||
response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px)
|
||
response = to_cpu(response)
|
||
except Exception:
|
||
response = -gaussian_laplace(filled, sigma=sigma_px)
|
||
else:
|
||
response = -gaussian_laplace(filled, sigma=sigma_px)
|
||
response[nan_mask] = np.nan
|
||
valid = response[~nan_mask]
|
||
std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
|
||
response = response / std_val
|
||
wavelet_stack.append(response)
|
||
|
||
# Gabor (directional) — 4 orientations at 3 key scales
|
||
gabor_scales_m = [5, 10, 20] # Key archaeological scales for linear features
|
||
gabor_orientations = [0, np.pi/4, np.pi/2, 3*np.pi/4] # 0°, 45°, 90°, 135°
|
||
gabor_scale_weights = {5: 2.0, 10: 1.8, 20: 1.5}
|
||
|
||
for scale_m in gabor_scales_m:
|
||
sigma_px = max(2, scale_m / resolution / 3) # Sigma relative to wavelength
|
||
wavelength_px = max(3, scale_m / resolution)
|
||
kernel_size = max(5, int(wavelength_px * 2.5))
|
||
if kernel_size % 2 == 0:
|
||
kernel_size += 1
|
||
|
||
for theta in gabor_orientations:
|
||
kernel = _gabor_kernel_2d(kernel_size, sigma_px, wavelength_px, theta)
|
||
# Normalize kernel
|
||
kernel = kernel / max(np.abs(kernel).max(), 1e-10)
|
||
response = convolve(filled, kernel, mode='constant', cval=0)
|
||
response[nan_mask] = np.nan
|
||
valid = response[~nan_mask]
|
||
std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
|
||
response = response / std_val
|
||
wavelet_stack.append(response)
|
||
|
||
# Build weights: Mexican Hat weights + Gabor weights
|
||
gabor_weights_list = []
|
||
for scale_m in gabor_scales_m:
|
||
w = gabor_scale_weights.get(scale_m, 1.0)
|
||
gabor_weights_list.extend([w] * len(gabor_orientations))
|
||
gabor_weights = np.array(gabor_weights_list)
|
||
all_weights = np.concatenate([mex_weights, gabor_weights])
|
||
|
||
# Weighted RMS combination
|
||
stack = np.array(wavelet_stack)
|
||
weights_3d = all_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(all_weights))
|
||
combined[nan_mask] = np.nan
|
||
|
||
_save_tif(output, combined.astype(np.float32), transform, crs)
|
||
logger.info(f" ✓ Ondelette terminée ({time.time()-t0:.1f}s){' [GPU]' 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
|
||
# ============================================================
|
||
|
||
def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Flow Accumulation — priority-flood sink filling + D8 accumulation (GPU if available).
|
||
|
||
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.
|
||
|
||
Uses log10 transformation for visualization.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Accumulation d'écoulement (flow accumulation){gpu_tag}...")
|
||
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)
|
||
|
||
# Priority-flood sink filling (Wang & Liu 2006, O(n log n))
|
||
from heapq import heappush, heappop
|
||
rows, cols = filled.shape
|
||
dem_filled = filled.copy()
|
||
|
||
# Use heap for priority-flood
|
||
visited = np.zeros((rows, cols), dtype=bool)
|
||
heap = []
|
||
# Seed with all border cells
|
||
for x in range(cols):
|
||
heappush(heap, (dem_filled[0, x], 0, x))
|
||
heappush(heap, (dem_filled[rows-1, x], rows-1, x))
|
||
visited[0, x] = True
|
||
visited[rows-1, x] = True
|
||
for y in range(1, rows-1):
|
||
heappush(heap, (dem_filled[y, 0], y, 0))
|
||
heappush(heap, (dem_filled[y, cols-1], y, cols-1))
|
||
visited[y, 0] = True
|
||
visited[y, cols-1] = True
|
||
|
||
while heap:
|
||
elev, cy, cx = heappop(heap)
|
||
for dy, dx in [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]:
|
||
ny, nx = cy + dy, cx + dx
|
||
if 0 <= ny < rows and 0 <= nx < cols and not visited[ny, nx]:
|
||
if dem_filled[ny, nx] > elev:
|
||
dem_filled[ny, nx] = elev
|
||
heappush(heap, (dem_filled[ny, nx], ny, nx))
|
||
visited[ny, nx] = True
|
||
|
||
logger.info(f" ✓ Sink filling terminé ({time.time()-t0:.1f}s)")
|
||
|
||
# D8 flow direction and accumulation
|
||
flow_acc = np.zeros((rows, cols), dtype=np.int64)
|
||
flow_dir = np.zeros((rows, cols), dtype=np.int8) - 1
|
||
|
||
# D8 neighbors (ordered by angle)
|
||
neighbors = [
|
||
(-1, 0), (-1, 1), ( 0, 1), ( 1, 1),
|
||
( 1, 0), ( 1, -1), ( 0, -1), (-1, -1)
|
||
]
|
||
|
||
# Process cells in ascending elevation order for correct accumulation
|
||
flat_idx = np.argsort(dem_filled.ravel())
|
||
flow_acc_flat = flow_acc.ravel()
|
||
|
||
for idx in flat_idx:
|
||
y = idx // cols
|
||
x = idx % cols
|
||
|
||
# Find steepest downslope neighbor
|
||
max_slope = -np.inf
|
||
best_dir = -1
|
||
for di, (dy, dx) in enumerate(neighbors):
|
||
ny, nx = y + dy, x + dx
|
||
if 0 <= ny < rows and 0 <= nx < cols:
|
||
slope = (dem_filled[y, x] - dem_filled[ny, nx])
|
||
dist = math.sqrt(dy**2 + dx**2)
|
||
slope_per_m = slope / dist
|
||
if slope_per_m > max_slope:
|
||
max_slope = slope_per_m
|
||
best_dir = di
|
||
|
||
flow_dir[y, x] = best_dir
|
||
flow_acc[y, x] = 1 # Count self
|
||
|
||
# Accumulate flow (upstream to downstream)
|
||
# Process in reverse elevation order (highest first)
|
||
for idx in reversed(flat_idx):
|
||
y = idx // cols
|
||
x = idx % cols
|
||
d = flow_dir[y, x]
|
||
if d >= 0:
|
||
dy, dx = neighbors[d]
|
||
ny, nx = y + dy, x + dx
|
||
if 0 <= ny < rows and 0 <= nx < cols:
|
||
flow_acc[ny, nx] += flow_acc[y, x]
|
||
|
||
logger.info(f" ✓ D8 accumulation terminé ({time.time()-t0:.1f}s)")
|
||
|
||
# Log transform for visualization
|
||
flow_result = np.log10(np.maximum(flow_acc.astype(np.float32), 1.0))
|
||
flow_result[nan_mask] = np.nan
|
||
|
||
_save_tif(output, flow_result, transform, crs)
|
||
logger.info(f" ✓ Flow accumulation terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur flow accumulation: {e}", exc_info=True)
|
||
return None
|
||
|
||
|
||
# ============================================================
|
||
# Anisotropic Openness
|
||
# ============================================================
|
||
|
||
def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None):
|
||
"""Anisotropic Openness - weighted directional openness, multi-radius (GPU if available).
|
||
|
||
Computes positive and negative openness with anisotropic weighting:
|
||
NW/SE directions weighted more heavily to enhance detection of structures
|
||
aligned NE-SW (common in French archaeological sites: villas, enclosures).
|
||
|
||
Multi-radius (25, 50, 100m) with equal weight, results std-normalized
|
||
for cross-tile comparability.
|
||
"""
|
||
gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
|
||
logger.info(f" → Openness Anisotropique (multi-rayon){gpu_tag}...")
|
||
t0 = time.time()
|
||
output = vis_dir / f"{basename}_aniso_open.tif"
|
||
|
||
try:
|
||
dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
|
||
_prepare_dem_for_raycast(dem_file, shared, resolution)
|
||
|
||
n_dirs = 8
|
||
# Anisotropic weights: emphasize NW-SE and NE-SW directions
|
||
weights = np.array([1.0, 1.5, 1.0, 1.5, 1.0, 1.5, 1.0, 1.5])
|
||
|
||
radii_m = [25, 50, 100]
|
||
max_dist = min(int(100 / res), 300)
|
||
|
||
pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m)
|
||
|
||
# Weighted combination across directions and radii
|
||
weight_total = np.sum(weights)
|
||
n_radii = len(radii_m)
|
||
|
||
pos_combined = np.zeros((rows, cols), dtype=np.float64)
|
||
neg_combined = np.zeros((rows, cols), dtype=np.float64)
|
||
|
||
for r_idx in range(n_radii):
|
||
for d_idx in range(n_dirs):
|
||
w = weights[d_idx]
|
||
pos_combined += pos_angles[d_idx, r_idx] * w / (n_radii * weight_total)
|
||
neg_combined += neg_angles[d_idx, r_idx] * w / (n_radii * weight_total)
|
||
|
||
aniso_result = np.degrees(pos_combined - neg_combined).astype(np.float32)
|
||
aniso_result[nan_mask] = np.nan
|
||
|
||
# Std normalization for cross-tile comparability
|
||
valid = aniso_result[~nan_mask]
|
||
if len(valid) > 0:
|
||
std_val = max(np.nanstd(valid), 0.01)
|
||
aniso_result = aniso_result / std_val
|
||
|
||
_save_tif(output, aniso_result, transform, crs)
|
||
logger.info(f" ✓ Openness anisotropique terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
|
||
return output
|
||
except Exception as e:
|
||
logger.error(f" ✗ Erreur openness anisotropique: {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, normalized to z-score, and contributes to the composite
|
||
layer_configs = [
|
||
# (filename_pattern, weight)
|
||
("mslrm", 2.5), # Multi-scale relief — strongest signal
|
||
("negative_openness", 1.8), # Fossés, dolines
|
||
("roughness", 1.5), # Surface irregularity
|
||
("wavelet", 1.5), # Circular + linear structures
|
||
("svf", 1.3), # Sky-view depressions
|
||
("flow_acc", 1.2), # Drainage channels / ditches
|
||
("aniso_open", 1.0), # Anisotropic structures
|
||
("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)
|
||
# Z-score normalization
|
||
valid = data[~nan_mask]
|
||
if len(valid) == 0:
|
||
continue
|
||
mean_val = np.nanmean(valid)
|
||
std_val = max(np.nanstd(valid), 0.01)
|
||
zscore = (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 = 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 RMS combination (unsigned — all deviations are suspicious)
|
||
combined = np.zeros((rows, cols), dtype=np.float64)
|
||
total_weight = 0.0
|
||
for layer_data, weight in layers:
|
||
combined += (layer_data ** 2) * weight
|
||
total_weight += weight
|
||
|
||
if total_weight > 0:
|
||
combined = np.sqrt(combined / total_weight)
|
||
|
||
# Apply n_sigma threshold: pixels below n_sigma are suppressed
|
||
combined = np.maximum(combined - n_sigma, 0.0)
|
||
|
||
# Rescale to 0–1 for visualization (percentile-based stretch)
|
||
valid_combined = combined[~nan_mask]
|
||
if len(valid_combined) > 0 and np.nanmax(valid_combined) > 0:
|
||
p99 = np.percentile(valid_combined, 99)
|
||
if p99 > 0:
|
||
combined = np.clip(combined / p99, 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
|