Refine MSRM: add small scales (3,8,15m), drop 100/200m, clip |z| to 3.0

Small archaeological features (ditches, walls, post-holes) were drowned out
by large-scale topography (100-200m). Fix: replace 100/200m with finer
scales (3, 8, 15m), weight 5-10m heaviest, clip absolute z-score to 3.0
before combination to prevent large-scale outliers from dominating.
This commit is contained in:
Antoine Jacquin
2026-06-01 23:18:15 +02:00
parent 8478106e51
commit 4dafae6d02

View File

@ -665,12 +665,16 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
# Adaptive scales: finer at higher resolution # Adaptive scales: finer at higher resolution
min_scale = max(2.0, resolution * 4) min_scale = max(2.0, resolution * 4)
candidate_scales = [2, 5, 10, 20, 50, 100, 200] # 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] sigmas = [s for s in candidate_scales if s >= min_scale]
# Archaeological weights: favor 5-25m range (ditches, enclosures, tumulus) # Weights: favor small-to-medium scales where archaeo features live
scale_weights = { scale_weights = {
2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6, 200: 0.4, 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]) weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
@ -685,10 +689,12 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px) local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
lrm = dem_np - local_mean lrm = dem_np - local_mean
lrm[nan_mask] = np.nan lrm[nan_mask] = np.nan
# Std normalization: x / std — preserves sign and contrast better than z-score # Std normalization
valid_lrm = lrm[~nan_mask] valid_lrm = lrm[~nan_mask]
lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
lrm = lrm / lrm_std 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)) lrm_stack.append(lrm.astype(np.float32))
# Weighted combination — preserve sign for RdBu_r colormap # Weighted combination — preserve sign for RdBu_r colormap