Remove LRM, TPI, aspect, curvature, paths + add flow accumulation, directional Gabor wavelets, multi-radius ray-tracing

lidar_pipeline/pipeline.py

@@ -57,11 +57,12 @@ _file_filter = FilePrefixFilter()
 from .dtm import classify_ground, create_dtm_fast
 from .visualizations import (
     SharedDEM,
-    generate_hillshade, generate_slope, generate_aspect, generate_curvature,
+    generate_hillshade, generate_slope,
-    generate_lrm, generate_openness,
+    generate_openness,
-    generate_mslrm, generate_tpi, generate_sailore,
+    generate_mslrm, generate_sailore,
     generate_roughness, generate_wavelet,
-    generate_svf, generate_aniso_open, generate_paths,
+    generate_svf, generate_aniso_open,
+    generate_flow_accumulation,
 )
 from .gpu import gpu_cleanup, num_gpus, restrict_gpus, safe_gpu_call
 from .ign import generate_ign_overlay
@@ -74,19 +75,15 @@ from .rendering import tif_to_png
 VIZ_STEPS = [
     ('hillshade',  generate_hillshade),
     ('slope',      generate_slope),
-    ('aspect',     generate_aspect),
+    ('mslrm',      generate_mslrm),
-    ('curvature',  generate_curvature),
+    ('sailore',    generate_sailore),
-    ('lrm',        generate_lrm),
     ('pos_open',   lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=True, shared=shared)),
     ('neg_open',   lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=False, shared=shared)),
-    ('mslrm',      generate_mslrm),
-    ('tpi',        generate_tpi),
-    ('sailore',    generate_sailore),
-    ('roughness',  generate_roughness),
     ('svf',        generate_svf),
     ('aniso_open', generate_aniso_open),
-    ('paths',      generate_paths),
+    ('roughness',  generate_roughness),
     ('wavelet',    generate_wavelet),
+    ('flow_acc',   generate_flow_accumulation),
     ('ortho',      lambda d, b, v, r: generate_ign_overlay(
                     d, b, v, r,
                     layer='ORTHOIMAGERY.ORTHOPHOTOS',

lidar_pipeline/rendering.py

@@ -81,13 +81,6 @@ _FRANCE_OUTLINE_L93 = np.array([
 COLORMAPS = {
     # === Famille RELIEF : rouge=surélévation, bleu=dépression ===
     # Diverging: rouge vif=positif, bleu vif=négatif, blanc=plat
-    'curvature': {
-        'cmap': 'bwr',
-        'title': 'Courbure (Convexité/Concavité du terrain)',
-        'legend': 'Changement de pente (1/m)\nRouge = Convexe (sommet de mur, levée)\nBleu = Concave (fond de fossé, dépression)',
-        'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées',
-        'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
-    },
     'mslrm': {
         'cmap': 'seismic',
         'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)',
@@ -95,20 +88,6 @@ COLORMAPS = {
         'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
-    'lrm': {
-        'cmap': 'seismic',
-        'title': 'LRM - Local Relief Model (échelle unique 15m)',
-        'legend': 'Écart local par rapport au terrain moyen (m)\nRouge = Surélévation (+{vmax:.2f}m)\nBleu = Dépression ({vmin:.2f}m)\nNoyau gaussien unique de 15m',
-        'description': 'Micro-relief à 15m seulement — voir MSRM pour toutes les échelles',
-        'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
-    },
-    'tpi': {
-        'cmap': 'seismic',
-        'title': 'TPI - Topographic Position Index (4 échelles)',
-        'legend': 'Position dans le paysage\nRouge = Plus haut que le voisinage (crête, plateau)\nBleu = Plus bas que le voisinage (fossé, vallée)\nCombine 4 échelles : 3m, 15m, 50m, 200m',
-        'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle',
-        'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
-    },
     'sailore': {
         'cmap': 'seismic',
         'title': 'SAILORE - LRM Auto-Adaptatif',
@@ -123,14 +102,6 @@ COLORMAPS = {
         'description': 'Openness avec pondération anisotropique — détecte mieux les structures alignées NW-SE et NE-SW',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
-    'paths': {
-        'cmap': 'inferno',
-        'title': 'Cheminement (chemins et sentiers)',
-        'legend': 'Openness directionnelle maximale\nJaune/vif = Chemin ou sentier\nNoir = Terrain plat\n\nDétecte les structures linéaires dans toutes les directions',
-        'description': 'Différence openness positive-négative maximale sur 8 directions — chemins, sentiers, ornières',
-        'vmin_mode': 'percentile', 'vmin_pct': 5,
-        'vmax_mode': 'percentile', 'vmax_pct': 99,
-    },
     # === Famille OUVERTURE : séquentiel, toujours positif ===
     'positive_openness': {
         'cmap': 'YlOrBr',
@@ -160,7 +131,7 @@ COLORMAPS = {
     'hillshade': {
         'cmap': 'gray',
         'title': 'Hillshade Multidirectionnel',
-        'legend': 'Illumination combinée de 6 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre',
+        'legend': 'Illumination combinée de 8 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre',
         'description': 'Ombres portées révélant micro-relief (murs, fossés, terrasses)',
         'vmin_mode': 'percentile', 'vmin_pct': 1,
         'vmax_mode': 'percentile', 'vmax_pct': 99,
@@ -173,14 +144,6 @@ COLORMAPS = {
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 97,
     },
-    'aspect': {
-        'cmap': 'twilight',
-        'title': 'Aspect (Direction des pentes)',
-        'legend': 'Direction vers laquelle le terrain descend\nCycle continu : Nord→Est→Sud→Ouest→Nord\nCouleurs perceptuellement uniformes (pas de saut de teinte)',
-        'description': 'Orientation des pentes — utile pour distinguer structures des formes naturelles',
-        'vmin_mode': 'fixed', 'vmin_val': 0,
-        'vmax_mode': 'fixed', 'vmax_val': 360,
-    },
     'roughness': {
         'cmap': 'plasma',
         'title': 'Rugosité Multi-Échelle (3m + 15m)',
@@ -191,11 +154,19 @@ COLORMAPS = {
     },
     'wavelet': {
         'cmap': 'cividis',
-        'title': 'Ondelette Mexican Hat (CWT multi-échelle)',
+        'title': 'Ondelette Mexican Hat + Gabor directionnelle (multi-échelle)',
-        'legend': 'Réponse de la transformée en ondelette\nÉchelles adaptées à la résolution\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
+        'legend': 'Réponse combinée Mexican Hat (circulaire) + Gabor (linéaire)\nÉchelles adaptées à la résolution\n4 orientations Gabor : 0°, 45°, 90°, 135°\n\nClair = Structure détectée\nSombre = Pas de structure',
-        'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires',
+        'description': 'Mexican Hat pour tumulus/enclos + Gabor pour chemins/murs/fossés',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
+    'flow_acc': {
+        'cmap': 'YlGn',
+        'title': 'Accumulation d\'Écoulement (Flow Accumulation)',
+        'legend': 'Log10 du nombre de cellules amont\nJaune = Fort accumulation (fossé, chenal, drainage)\nVert clair = Faible accumulation\n\nDétection fossés et linéaires hydrologiques',
+        'description': 'Priority-flood + D8 — détecte fossés archéologiques et drainages',
+        'vmin_mode': 'fixed', 'vmin_val': 0,
+        'vmax_mode': 'percentile', 'vmax_pct': 98,
+    },
 }
 
 # RGB entries (ortho/topo) are handled specially
@@ -842,8 +813,7 @@ def generate_pdf_report(basename, vis_dir, pdf_dir, resolution):
     # Sort analysis files by archaeological priority
     order = ['mslrm', 'svf', 'negative_openness',
               'positive_openness', 'aniso_open', 'sailore', 'hillshade_multi',
-              'lrm', 'tpi', 'slope', 'curvature', 'aspect',
+              'flow_acc', 'slope', 'roughness', 'wavelet']
-              'roughness', 'wavelet']
 
     def sort_key(f):
         name = f.stem.lower()

lidar_pipeline/tests/test_visualizations.py

@@ -47,52 +47,8 @@ class TestSlope:
             assert np.nanmax(data) <= 90
 
 
-class TestAspect:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_aspect
-        result = generate_aspect(synthetic_dem, "test", tmp_output_dir, 5.0)
-        assert result is not None
-        assert result.exists()
-
-    def test_aspect_values_0_360(self, synthetic_dem, tmp_output_dir):
-        import rasterio
-        from lidar_pipeline.visualizations import generate_aspect
-        result = generate_aspect(synthetic_dem, "test", tmp_output_dir, 5.0)
-        with rasterio.open(result) as src:
-            data = src.read(1)
-            valid = data[~np.isnan(data)]
-            assert np.nanmin(valid) >= 0
-            assert np.nanmax(valid) <= 360
-
-
-class TestCurvature:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_curvature
-        result = generate_curvature(synthetic_dem, "test", tmp_output_dir, 5.0)
-        assert result is not None
-        assert result.exists()
-
-
 # --- GPU-accelerated visualizations ---
 
-class TestLRM:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_lrm
-        result = generate_lrm(synthetic_dem, "test", tmp_output_dir, 5.0)
-        assert result is not None
-        assert result.exists()
-
-    def test_lrm_has_positive_negative(self, synthetic_dem, tmp_output_dir):
-        import rasterio
-        from lidar_pipeline.visualizations import generate_lrm
-        result = generate_lrm(synthetic_dem, "test", tmp_output_dir, 5.0)
-        with rasterio.open(result) as src:
-            data = src.read(1)
-            # LRM should have both positive and negative values
-            assert np.nanmax(data) > 0
-            assert np.nanmin(data) < 0
-
-
 class TestSVF:
     def test_generates_tif(self, synthetic_dem, tmp_output_dir):
         from lidar_pipeline.visualizations import generate_svf
@@ -133,15 +89,6 @@ class TestMSLRM:
         assert result.exists()
 
 
-class TestTPI:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_tpi
-        result = generate_tpi(synthetic_dem, "test", tmp_output_dir, 5.0)
-        assert result is not None
-        assert result.exists()
-
-
-
 class TestSAILORE:
     def test_generates_tif(self, synthetic_dem, tmp_output_dir):
         from lidar_pipeline.visualizations import generate_sailore
@@ -167,14 +114,6 @@ class TestRoughness:
             assert np.nanmin(data) >= 0
 
 
-class TestAnomalies:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_anomalies
-        result = generate_anomalies(synthetic_dem, "test", tmp_output_dir, 5.0)
-        assert result is not None
-        assert result.exists()
-
-
 class TestWavelet:
     def test_generates_tif(self, synthetic_dem, tmp_output_dir):
         from lidar_pipeline.visualizations import generate_wavelet
@@ -183,50 +122,38 @@ class TestWavelet:
         assert result.exists()
 
 
-class TestFlow:
+class TestFlowAccumulation:
     def test_generates_tif(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_flow
+        from lidar_pipeline.visualizations import generate_flow_accumulation
-        result = generate_flow(synthetic_dem, "test", tmp_output_dir, 5.0)
+        result = generate_flow_accumulation(synthetic_dem, "test", tmp_output_dir, 5.0)
         assert result is not None
         assert result.exists()
 
     def test_flow_log_values(self, synthetic_dem, tmp_output_dir):
         import rasterio
-        from lidar_pipeline.visualizations import generate_flow
+        from lidar_pipeline.visualizations import generate_flow_accumulation
-        result = generate_flow(synthetic_dem, "test", tmp_output_dir, 5.0)
+        result = generate_flow_accumulation(synthetic_dem, "test", tmp_output_dir, 5.0)
         with rasterio.open(result) as src:
             data = src.read(1)
-            # log1p(x) >= 0 for x >= 0
+            # log10(x) >= 0 for x >= 1
             valid = data[~np.isnan(data)]
             assert np.nanmin(valid) >= 0
 
 
-class TestLocalDominance:
+class TestRayTrace:
-    def test_generates_tif(self, synthetic_dem, tmp_output_dir):
+    def test_rays_are_traced(self, synthetic_dem, tmp_output_dir):
-        from lidar_pipeline.visualizations import generate_local_dominance
+        """Verify _ray_trace_horizons returns expected shapes."""
-        result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0)
+        from lidar_pipeline.visualizations import _ray_trace_horizons, _prepare_dem_for_raycast
-        assert result is not None
-        assert result.exists()
-        assert result.suffix == ".tif"
-
-    def test_dominance_values_0_1(self, synthetic_dem, tmp_output_dir):
         import rasterio
-        from lidar_pipeline.visualizations import generate_local_dominance
+        with rasterio.open(synthetic_dem) as src:
-        result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0)
+            dem_np = src.read(1)
-        with rasterio.open(result) as src:
+        rows, cols = dem_np.shape
-            data = src.read(1)
+        # Create a simple filled DEM for testing
-            valid = data[~np.isnan(data)]
+        import numpy as np
-            assert np.nanmin(valid) >= 0, "Local dominance should be >= 0"
+        filled = np.nan_to_num(dem_np, nan=0)
-            assert np.nanmax(valid) <= 1, "Local dominance should be <= 1"
+        # Test with numpy (no GPU)
-
+        pos, neg = _ray_trace_horizons(
-    def test_dominance_nan_mask_preserved(self, synthetic_dem, tmp_output_dir):
+            filled, rows, cols, 5.0, n_dirs=4, max_dist=10, radii_m=[25, 50]
-        """Check that NaN zones from original DEM are preserved."""
+        )
-        import rasterio
+        assert pos.shape == (4, 2, rows, cols)
-        from lidar_pipeline.visualizations import generate_local_dominance
+        assert neg.shape == (4, 2, rows, cols)
-        result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0)
-        # The synthetic DEM has no NaN, so this just verifies the output is valid
-        with rasterio.open(result) as src:
-            data = src.read(1)
-            # Shape should match input
-            assert data.shape[0] > 0
-            assert data.shape[1] > 0