Fix corrupted COPC detection, add CSF→SMRF fallback, improve MSRM colormap, add SVF and anisotropic openness

- validate_laz: verify point data accessibility (not just headers) to detect
  corrupted COPC files that pass header checks but fail on data reads
- classify_ground: fallback from CSF to SMRF when CSF produces no ground
  points or PDAL errors (fixes 2/9 failing tiles)
- MSRM: preserve sign in weighted combination so RdBu_r colormap shows
  both red (elevated) and blue (depressed) instead of red only
- Add Sky-View Factor (SVF) visualization: cos²(horizon angle) over 16
  directions, excellent for archaeological earthwork detection
- Add Anisotropic Openness: directional weighting (NW-SE/NE-SW) enhances
  linear feature detection aligned with common settlement patterns
- Remove anomalies and flow visualizations (replaced by SVF + aniso_open)
- Location inset: use IGN topographic map at zoom 10 instead of simplified
  France outline, with red rectangle marker and fallback
- Remove flow (hydrological accumulation) from VIZ_STEPS

lidar_pipeline/dtm.py

@@ -125,15 +125,15 @@ def create_csf_pipeline(input_laz, output_las):
 
 
 def validate_laz(laz_file):
-    """Quick integrity check for a LAZ/LAS file.
+    """Integrity check for a LAZ/LAS file.
 
-    Tries laspy first (fast header read), then PDAL as fallback for COPC files
+    Verifies that both the header AND point data are readable. Some corrupted
-    that laspy cannot read. Also checks that the file contains points.
+    COPC files have valid headers but inaccessible point data (LazrsError:
+    failed to fill whole buffer). Such files must be re-downloaded.
 
     Returns:
-        True if file is readable and contains points, False otherwise.
+        True if file is readable and contains accessible points, False otherwise.
     """
-    # Try laspy first (fast)
     import laspy
     try:
         with laspy.open(str(laz_file)) as f:
@@ -143,6 +143,15 @@ def validate_laz(laz_file):
                 logger.error(f"  ✗ Fichier vide (0 points): {laz_file.name}")
                 logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
                 return False
+            # Verify point data is actually accessible (not just header metadata)
+            try:
+                for _ in f.chunk_iterator(1000):
+                    break  # Read just one chunk to confirm data integrity
+            except Exception as chunk_err:
+                logger.error(f"  ✗ Données inaccessibles (fichier corrompu?): {laz_file.name}")
+                logger.error(f"    Erreur: {chunk_err}")
+                logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
+                return False
         return True
     except Exception:
         pass
@@ -166,6 +175,18 @@ def validate_laz(laz_file):
                     return False
             except Exception:
                 pass  # Can't parse — assume valid
+            # Verify PDAL can actually read point data (not just header)
+            try:
+                test_result = subprocess.run(
+                    ["pdal", "info", str(laz_file), "--point", "1"],
+                    capture_output=True, text=True, timeout=60
+                )
+                if test_result.returncode != 0:
+                    logger.error(f"  ✗ Données inaccessibles (PDAL): {laz_file.name}")
+                    logger.error(f"    → Re-télécharger depuis https://ign.fr/lidar-hd")
+                    return False
+            except (subprocess.TimeoutExpired, FileNotFoundError):
+                pass  # Timeout — assume valid, will fail later if corrupted
             return True
         logger.error(f"  ✗ Fichier illisible: {laz_file.name}")
         logger.error(f"    PDAL: {result.stderr.strip()[:200]}")
@@ -347,6 +368,9 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
         if output_las.exists() and output_las.stat().st_size < 100:
             logger.error(f"  ✗ Fichier ground vide (taille < 100 octets)")
             output_las.unlink(missing_ok=True)
+            # Fallback: if CSF produced no ground points, retry with SMRF
+            if method == 'csf':
+                return _fallback_to_smrf(laz_file, temp_dir, laz_base, force)
             return None
         logger.info(f"  ✓ Classification sol {method.upper()} terminée")
         return output_las
@@ -354,6 +378,10 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
         error_msg = e.stderr.decode() if e.stderr else str(e)
         logger.warning(f"  ✗ Erreur classification PDAL ({method.upper()}): {error_msg}")
 
+        # Fallback: if CSF failed, retry with SMRF
+        if method == 'csf':
+            return _fallback_to_smrf(laz_file, temp_dir, laz_base, force)
+
         # Try repairing file with laspy if PDAL fails on EVLR/VLR
         if 'VLR' in error_msg or 'Invalid' in error_msg:
             logger.info(f"  → Tentative de réparation du fichier avec laspy...")
@@ -378,6 +406,58 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
         return None
 
 
+def _fallback_to_smrf(laz_file, temp_dir, laz_base, force=False):
+    """Retry ground classification with SMRF when CSF fails.
+
+    CSF (Cloth Simulation Filter) can fail on certain terrain types where
+    SMRF (Simple Morphological Filter) succeeds. This fallback ensures
+    processing continues even when auto-detection selects CSF incorrectly.
+
+    Args:
+        laz_file: Path to input LAZ/LAS file.
+        temp_dir: Directory for temporary files.
+        laz_base: Base name for the file.
+        force: If True, reclassify even if output exists.
+
+    Returns:
+        Path to classified ground LAS file, or None on failure.
+    """
+    logger.info(f"  → Basculement CSF → SMRF (fallback)")
+
+    # Clean up failed CSF output if it exists
+    csf_output = temp_dir / f"{laz_base}_ground_csf.las"
+    if csf_output.exists():
+        csf_output.unlink(missing_ok=True)
+
+    output_las = temp_dir / f"{laz_base}_ground_smrf.las"
+
+    if output_las.exists() and not force:
+        logger.info(f"  Classification SMRF déjà existante — fichier réutilisé")
+        return output_las
+
+    pipeline_json = _create_ground_pipeline(laz_file, output_las, 'smrf')
+    pipeline_file = temp_dir / "pipeline_smrf.json"
+
+    with open(pipeline_file, 'w') as f:
+        f.write(pipeline_json)
+
+    try:
+        subprocess.run(
+            ["pdal", "pipeline", str(pipeline_file)],
+            capture_output=True, check=True
+        )
+        if output_las.exists() and output_las.stat().st_size < 100:
+            logger.error(f"  ✗ Fichier ground SMRF vide (taille < 100 octets)")
+            output_las.unlink(missing_ok=True)
+            return None
+        logger.info(f"  ✓ Classification sol SMRF terminée (fallback depuis CSF)")
+        return output_las
+    except subprocess.CalledProcessError as e:
+        error_msg = e.stderr.decode() if e.stderr else str(e)
+        logger.error(f"  ✗ Échec classification SMRF (fallback): {error_msg}")
+        return None
+
+
 def _repair_laz_with_laspy(input_laz, output_las):
     """Try to repair a corrupt LAZ file by re-reading with laspy and saving as LAS.
 

lidar_pipeline/pipeline.py

@@ -60,8 +60,8 @@ from .visualizations import (
     generate_hillshade, generate_slope, generate_aspect, generate_curvature,
     generate_lrm, generate_openness,
     generate_mslrm, generate_tpi, generate_sailore,
-    generate_roughness, generate_anomalies, generate_wavelet,
+    generate_roughness, generate_wavelet,
-    generate_flow,
+    generate_svf, generate_aniso_open,
 )
 from .gpu import gpu_cleanup
 from .ign import generate_ign_overlay
@@ -83,9 +83,9 @@ VIZ_STEPS = [
     ('tpi',        generate_tpi),
     ('sailore',    generate_sailore),
     ('roughness',  generate_roughness),
-    ('anomalies',  generate_anomalies),
+    ('svf',        generate_svf),
+    ('aniso_open', generate_aniso_open),
     ('wavelet',    generate_wavelet),
-    ('flow',       generate_flow),
     ('ortho',      lambda d, b, v, r: generate_ign_overlay(
                     d, b, v, r,
                     layer='ORTHOIMAGERY.ORTHOPHOTOS',

lidar_pipeline/rendering.py

@@ -21,12 +21,14 @@ try:
 except ImportError:
     HAS_WARP = False
 
+# Cache for IGN location map tiles (avoid re-downloading for each visualization)
+_location_map_cache = {}
+
 import matplotlib
 matplotlib.use('Agg')
 import matplotlib.pyplot as plt
 from matplotlib import rcParams
-from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch
+from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch, FancyBboxPatch
-from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 
 rcParams['figure.dpi'] = 150
 rcParams['savefig.dpi'] = 300
@@ -35,6 +37,32 @@ rcParams['font.size'] = 10
 logger = logging.getLogger("lidar")
 
 
+# ============================================================
+# Simplified France outline in Lambert 93 (EPSG:2154)
+# Used for location inset map on each visualization
+# ============================================================
+_FRANCE_OUTLINE_L93 = np.array([
+    [109000, 6385000], [134000, 6410000], [153000, 6430000], [173000, 6445000],
+    [200000, 6460000], [250000, 6475000], [300000, 6490000], [350000, 6500000],
+    [400000, 6505000], [450000, 6510000], [500000, 6510000], [550000, 6510000],
+    [600000, 6505000], [650000, 6500000], [700000, 6495000], [750000, 6485000],
+    [800000, 6470000], [840000, 6460000], [880000, 6450000], [920000, 6435000],
+    [950000, 6425000], [980000, 6415000], [1010000, 6405000], [1040000, 6395000],
+    [1060000, 6385000], [1080000, 6370000], [1100000, 6355000], [1120000, 6340000],
+    [1140000, 6320000], [1160000, 6300000], [1175000, 6280000], [1185000, 6260000],
+    [1190000, 6240000], [1195000, 6220000], [1198000, 6200000], [1196000, 6180000],
+    [1192000, 6160000], [1185000, 6140000], [1175000, 6120000], [1160000, 6100000],
+    [1140000, 6085000], [1120000, 6070000], [1095000, 6060000], [1070000, 6050000],
+    [1040000, 6040000], [1000000, 6035000], [950000, 6035000], [900000, 6035000],
+    [850000, 6040000], [800000, 6045000], [750000, 6050000], [700000, 6055000],
+    [650000, 6060000], [600000, 6065000], [550000, 6070000], [500000, 6075000],
+    [450000, 6080000], [400000, 6085000], [350000, 6095000], [300000, 6110000],
+    [250000, 6125000], [200000, 6145000], [160000, 6170000], [130000, 6200000],
+    [110000, 6230000], [100000, 6260000], [95000, 6290000], [100000, 6310000],
+    [105000, 6340000], [109000, 6385000],
+])
+
+
 # ============================================================
 # Colormap registry
 # ============================================================
@@ -126,12 +154,20 @@ COLORMAPS = {
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 97,
     },
-    'anomalies': {
+    'svf': {
-        'cmap': 'coolwarm',
+        'cmap': 'gray_r',
-        'title': 'Anomalies Statistiques (MSRM multi-échelle + Moran\'s I)',
+        'title': 'Sky-View Factor (fraction de ciel visible)',
-        'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine MSRM normalisé (intensité) et\nMoran\'s I (regroupement spatial)',
+        'legend': 'Proportion de ciel visible depuis chaque point\nBlanc = Ciel dégagé (sommet, plateau, levée)\nNoir = Ciel masqué (vallée, fossé, tranchée)\nMoyenne de cos²(angle horizon) sur 16 directions',
-        'description': 'Détecte uniquement les anomalies statistiquement significatives — filtre le bruit de fond',
+        'description': 'Détection de micro-relief — fossés sombres, levées claires, complémentaire de l\'openness',
-        'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
+        'vmin_mode': 'fixed', 'vmin_val': 0,
+        'vmax_mode': 'fixed', 'vmax_val': 1,
+    },
+    'aniso_open': {
+        'cmap': 'RdBu_r',
+        'title': 'Openness Anisotropique (pondération directionnelle)',
+        'legend': 'Openness positive - négative pondérée (degrés)\nRouge = Surélévation dominante (mur, levée)\nBleu = Dépression dominante (fossé, doline)\nPondère les directions NW-SE et NE-SW davantage',
+        '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),
     },
     'wavelet': {
         'cmap': 'cividis',
@@ -231,6 +267,66 @@ def _apply_colormap(data, tif_file):
     return data, 'terrain', title, 'Altitude normalisée', '', False
 
 
+def _download_location_map(min_x, max_x, min_y, max_y):
+    """Download a wide-area IGN topographic map for location context.
+
+    Downloads a zoomed-out IGN PLANIGNV2 tile covering 5-10x the processed
+    zone extent, giving a wider geographic context. Results are cached to
+    avoid re-downloading for each visualization in the same tile.
+
+    Args:
+        min_x, max_x, min_y, max_y: DTM bounds in Lambert 93.
+
+    Returns:
+        numpy array (H, W, 3) uint8, or None on failure.
+    """
+    import hashlib
+
+    # Cache key based on rounded coordinates (1km grid)
+    cache_key = (round(min_x, -3), round(max_x, -3), round(min_y, -3), round(max_y, -3))
+    if cache_key in _location_map_cache:
+        return _location_map_cache[cache_key]
+
+    from .ign import download_ign_tiles, _optimal_zoom_level
+
+    if not HAS_WARP:
+        return None
+
+    try:
+        # Compute center coordinates for zoom calculation
+        center_x = (min_x + max_x) / 2
+        center_y = (min_y + max_y) / 2
+        clons, clats = warp_transform('EPSG:2154', 'EPSG:4326', [center_x], [center_y])
+        center_lat = clats[0]
+        center_lon = clons[0]
+
+        # Use a much lower zoom level for context (wider view)
+        # Zoom 10 gives ~150km per 256px tile — perfect for a small location map
+        context_zoom = 10
+
+        # Expand bounds by 3x in each direction for wider context
+        extent_x = max_x - min_x
+        extent_y = max_y - min_y
+        context_min_x = center_x - extent_x * 2
+        context_max_x = center_x + extent_x * 2
+        context_min_y = center_y - extent_y * 2
+        context_max_y = center_y + extent_y * 2
+
+        result = download_ign_tiles(
+            context_min_x, context_max_x, context_min_y, context_max_y,
+            layer='GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2',
+            zoom_level=context_zoom
+        )
+
+        if result is not None:
+            _location_map_cache[cache_key] = result
+
+        return result
+    except Exception as e:
+        logger.debug(f"    Carte de localisation IGN non disponible: {e}")
+        return None
+
+
 def _nice_scale(extent_m):
     """Choose a nice round scale distance that fits well in the image.
 
@@ -397,12 +493,16 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
 
             ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10)
 
-            # Colorbar/legend area — always at the same position for consistent layout
+            # Colorbar/legend area — reduced height to leave room for compass rose above
             cbar_left = data_left + data_width_frac + 0.02
             cbar_width = 0.04
+            compass_height = 0.07
+            compass_gap = 0.02
+            cbar_bottom = data_bottom
+            cbar_height = data_height_frac - compass_height - compass_gap
             if is_rgb:
                 # RGB: descriptive text label instead of gradient colorbar
-                cbar_ax = fig.add_axes([cbar_left, data_bottom, cbar_width, data_height_frac])
+                cbar_ax = fig.add_axes([cbar_left, cbar_bottom, cbar_width, cbar_height])
                 cbar_ax.set_xticks([])
                 cbar_ax.set_yticks([])
                 cbar_ax.text(0.5, 0.5, legend_label, transform=cbar_ax.transAxes,
@@ -411,7 +511,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
                              wrap=True)
                 cbar_ax.set_frame_on(False)
             elif is_rgba and saved_cmap is not None:
-                cbar_ax = fig.add_axes([cbar_left, data_bottom, cbar_width, data_height_frac])
+                cbar_ax = fig.add_axes([cbar_left, cbar_bottom, cbar_width, cbar_height])
                 sm = plt.cm.ScalarMappable(cmap=saved_cmap,
                                            norm=plt.Normalize(vmin=saved_vmin, vmax=saved_vmax))
                 sm.set_array([])
@@ -420,7 +520,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
                 cbar.outline.set_linewidth(1.5)
                 cbar.set_label(legend_label, fontsize=10, fontweight='bold')
             else:
-                cbar_ax = fig.add_axes([cbar_left, data_bottom, cbar_width, data_height_frac])
+                cbar_ax = fig.add_axes([cbar_left, cbar_bottom, cbar_width, cbar_height])
                 cbar = plt.colorbar(im, cax=cbar_ax)
                 cbar.ax.tick_params(labelsize=9, width=1.5)
                 cbar.outline.set_linewidth(1.5)
@@ -461,13 +561,16 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
                 spine.set_color('black')
                 spine.set_linewidth(0.8)
 
-            # North arrow — compass rose style
+            # North arrow — compass rose style, positioned above the colorbar
-            north_ax = inset_axes(ax, width="5%", height="9%", loc='upper right',
+            compass_bottom = data_bottom + data_height_frac + 0.02
-                                 bbox_to_anchor=(-0.03, 0.08, 1, 1), bbox_transform=ax.transAxes)
+            compass_height = 0.07
+            compass_width = cbar_width + 0.03
+            north_ax = fig.add_axes([cbar_left, compass_bottom, compass_width, compass_height])
             north_ax.set_xlim(-1.2, 1.2)
             north_ax.set_ylim(-0.5, 1.5)
             north_ax.axis('off')
             north_ax.set_aspect('equal')
+            north_ax.set_facecolor('white')
             # N arrow
             north_ax.annotate('N', xy=(0, 1.3), fontsize=11, fontweight='bold',
                              ha='center', va='bottom', color='#b22222')
@@ -566,6 +669,54 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
                         [bar_bottom_y - 0.05, bar_top_y + 0.05],
                         color='black', linewidth=1, transform=info_ax.transAxes, clip_on=False)
 
+            # Location inset map — IGN topographic background with processed zone marker
+            map_ax = fig.add_axes([0.84, 0.02, 0.14, 0.12])
+
+            # Try to download a wide-area IGN topo map for location context
+            location_map = _download_location_map(min_x, max_x, min_y, max_y)
+            if location_map is not None:
+                # Draw IGN topo map as background
+                map_ax.imshow(location_map, aspect='auto', extent=[
+                    min_x - (max_x - min_x) * 2, max_x + (max_x - min_x) * 2,
+                    min_y - (max_y - min_y) * 2, max_y + (max_y - min_y) * 2
+                ])
+                # Mark the processed zone with a red rectangle
+                rect_x1, rect_x2 = min_x, max_x
+                rect_y1, rect_y2 = min_y, max_y
+                map_ax.add_patch(RectPatch((rect_x1, rect_y1),
+                                           rect_x2 - rect_x1, rect_y2 - rect_y1,
+                                           facecolor='#ff3333', edgecolor='#cc0000',
+                                           linewidth=1.5, alpha=0.6, zorder=5))
+            else:
+                # Fallback: simplified France outline
+                map_ax.set_facecolor('#e8e8e8')
+                france = _FRANCE_OUTLINE_L93
+                map_ax.fill(france[:, 0] / 1000, france[:, 1] / 1000,
+                            facecolor='#f5f0e6', edgecolor='#888888', linewidth=0.8)
+                rect_x1, rect_x2 = min_x / 1000, max_x / 1000
+                rect_y1, rect_y2 = min_y / 1000, max_y / 1000
+                map_ax.add_patch(RectPatch((rect_x1, rect_y1),
+                                           rect_x2 - rect_x1, rect_y2 - rect_y1,
+                                           facecolor='#ff3333', edgecolor='#cc0000',
+                                           linewidth=1.2, alpha=0.7, zorder=5))
+                map_ax.set_xlim(france[:, 0].min() / 1000 - 50, france[:, 0].max() / 1000 + 50)
+                map_ax.set_ylim(france[:, 1].min() / 1000 - 50, france[:, 1].max() / 1000 + 50)
+
+            map_ax.set_aspect('equal')
+            map_ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
+            for spine in map_ax.spines.values():
+                spine.set_edgecolor('#aaaaaa')
+                spine.set_linewidth(0.5)
+            # Label with coordinates
+            if gps_coords:
+                nw_lat, nw_lon = gps_coords['NW']
+                se_lat, se_lon = gps_coords['SE']
+                map_ax.set_title(f"{nw_lat:.2f}°N  {nw_lon:.2f}°E",
+                                 fontsize=6, pad=1, color='#333333')
+            else:
+                map_ax.set_title(f"X:{min_x/1000:.0f} Y:{min_y/1000:.0f} km L93",
+                                 fontsize=6, pad=1, color='#333333')
+
             fig.patch.set_facecolor('white')
 
             # Save as PNG then convert to final format — fixed layout, no bbox_inches='tight'

lidar_pipeline/visualizations.py

@@ -702,13 +702,17 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
             lrm = lrm / lrm_std
             lrm_stack.append(lrm.astype(np.float32))
 
-        # Weighted combination
+        # 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')
-                mslrm = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
+                # 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_tag}")
@@ -970,8 +974,6 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
             valid_lrm = lrm[~nan_mask]
             lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
             lrm_norm = lrm / lrm_std
-            else:
-                lrm_norm = lrm
             lrm_stack.append(lrm_norm.astype(np.float32))
 
         # Weighted RMS combination (favor 5-25m scales)
@@ -1276,3 +1278,179 @@ def generate_flow(dem_file, basename, vis_dir, resolution, shared=None):
     except Exception as e:
         logger.error(f"  ✗ Erreur flux: {e}", exc_info=True)
         return None
+
+
+# ============================================================
+# Sky-View Factor (SVF)
+# ============================================================
+
+def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
+    """Sky-View Factor - fraction of sky visible from each point (GPU if available).
+
+    SVF = average of cos²(horizon_angle) across 8 directions.
+    High SVF (near 1) = open sky (ridgetop, plateau)
+    Low SVF (near 0) = enclosed sky (valley, deep trench)
+
+    Excellent for detecting archaeological earthworks: ditches appear dark,
+    embankments appear bright. Complements openness which uses raw angles.
+    """
+    gpu_tag = " [GPU]" if HAS_GPU else ""
+    logger.info(f"  → Sky-View Factor{gpu_tag}...")
+    t0 = time.time()
+    output = vis_dir / f"{basename}_svf.tif"
+
+    try:
+        if shared:
+            transform = shared.transform
+            crs = shared.crs
+            dem_np = shared.dem_np
+            rows, cols = dem_np.shape
+            res = resolution
+            dem = to_gpu(shared.filled) if HAS_GPU else shared.filled
+            nan_mask = shared.nan_mask
+        else:
+            dem_np, transform, crs = _read_dem(dem_file)
+            rows, cols = dem_np.shape
+            res = resolution
+            nan_mask = np.isnan(dem_np)
+            filled, _ = _fill_nans(dem_np)
+            dem = to_gpu(filled) if HAS_GPU else filled
+
+        n_dirs = 16  # More directions for smoother SVF
+        angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
+        dx_dir = np.cos(angles)
+        dy_dir = np.sin(angles)
+        max_dist = int(100 / res)
+
+        padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
+        svf_sum = xp.zeros_like(dem)
+
+        for d_idx in range(n_dirs):
+            ddx, ddy = dx_dir[d_idx], dy_dir[d_idx]
+            # Find maximum horizon elevation angle in this direction
+            max_horizon_angle = xp.zeros_like(dem)
+
+            for step in range(1, max_dist + 1):
+                px = int(round(ddx * step))
+                py = int(round(ddy * step))
+                dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
+                if dist_m < res * 0.5:
+                    continue
+
+                elev_diff = padded[max_dist + py:max_dist + py + rows,
+                                  max_dist + px:max_dist + px + cols] - dem
+
+                # Horizon angle from horizontal (positive = terrain above viewer)
+                angle = xp.arctan2(elev_diff, dist_m)
+                max_horizon_angle = xp.where(xp.isnan(angle), max_horizon_angle,
+                                             xp.maximum(max_horizon_angle, xp.nan_to_num(angle, nan=0)))
+
+            # SVF uses cos²(horizon angle) — fraction of visible sky in this direction
+            cos2 = xp.cos(max_horizon_angle) ** 2
+            svf_sum += cos2
+
+        svf_result = to_cpu(svf_sum / n_dirs).astype(np.float32)
+        svf_result[nan_mask] = np.nan
+        _save_tif(output, svf_result, transform, crs)
+        logger.info(f"  ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}")
+        return output
+    except Exception as e:
+        logger.error(f"  ✗ Erreur SVF: {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 emphasizing oblique directions (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).
+
+    The anisotropic weighting makes subtle linear features more visible than
+    standard isotropic openness which averages all directions equally.
+    """
+    gpu_tag = " [GPU]" if HAS_GPU else ""
+    logger.info(f"  → Openness Anisotropique{gpu_tag}...")
+    t0 = time.time()
+    output = vis_dir / f"{basename}_aniso_open.tif"
+
+    try:
+        if shared:
+            transform = shared.transform
+            crs = shared.crs
+            dem_np = shared.dem_np
+            rows, cols = dem_np.shape
+            res = resolution
+            dem = to_gpu(shared.filled) if HAS_GPU else shared.filled
+            nan_mask = shared.nan_mask
+        else:
+            dem_np, transform, crs = _read_dem(dem_file)
+            rows, cols = dem_np.shape
+            res = resolution
+            nan_mask = np.isnan(dem_np)
+            filled, _ = _fill_nans(dem_np)
+            dem = to_gpu(filled) if HAS_GPU else filled
+
+        n_dirs = 8
+        angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
+        dx_dir = np.cos(angles)
+        dy_dir = np.sin(angles)
+
+        # Anisotropic weights: emphasize NW-SE and NE-SW directions
+        # These orientations are most productive for detecting archaeological features
+        # aligned with Roman and medieval settlement patterns in France
+        weights = np.array([1.0, 1.5, 1.0, 1.5, 1.0, 1.5, 1.0, 1.5])
+
+        max_dist = int(100 / res)
+        padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
+
+        pos_sum = xp.zeros_like(dem)
+        neg_sum = xp.zeros_like(dem)
+        weight_total = 0.0
+
+        for d_idx in range(n_dirs):
+            ddx, ddy = dx_dir[d_idx], dy_dir[d_idx]
+            w = weights[d_idx]
+            weight_total += w
+
+            max_pos_angle = xp.zeros_like(dem)
+            max_neg_angle = xp.zeros_like(dem)
+
+            for step in range(1, max_dist + 1):
+                px = int(round(ddx * step))
+                py = int(round(ddy * step))
+                dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
+                if dist_m < res * 0.5:
+                    continue
+
+                elev_diff = padded[max_dist + py:max_dist + py + rows,
+                                  max_dist + px:max_dist + px + cols] - dem
+
+                # Positive openness: max zenith angle
+                pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m)
+                max_pos_angle = xp.where(xp.isnan(pos_angle), max_pos_angle,
+                                        xp.maximum(max_pos_angle, xp.nan_to_num(pos_angle, nan=0)))
+
+                # Negative openness: max nadir angle
+                neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m)
+                max_neg_angle = xp.where(xp.isnan(neg_angle), max_neg_angle,
+                                        xp.maximum(max_neg_angle, xp.nan_to_num(neg_angle, nan=0)))
+
+            pos_sum += max_pos_angle * w
+            neg_sum += max_neg_angle * w
+
+        # Combined: positive minus negative openness (anisotropic)
+        pos_avg = pos_sum / weight_total
+        neg_avg = neg_sum / weight_total
+        aniso_result = to_cpu(xp.degrees(pos_avg - neg_avg)).astype(np.float32)
+        aniso_result[nan_mask] = np.nan
+        _save_tif(output, aniso_result, transform, crs)
+        logger.info(f"  ✓ Openness anisotropique terminé ({time.time()-t0:.1f}s){gpu_tag}")
+        return output
+    except Exception as e:
+        logger.error(f"  ✗ Erreur openness anisotropique: {e}", exc_info=True)
+        return None