diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py index 9e7f15c..04b7a66 100644 --- a/lidar_pipeline/dtm.py +++ b/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 - that laspy cannot read. Also checks that the file contains points. + Verifies that both the header AND point data are readable. Some corrupted + 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. diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 0f4f23e..17814df 100644 --- a/lidar_pipeline/pipeline.py +++ b/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_flow, + generate_roughness, generate_wavelet, + 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', diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index ab6ebef..7c96aeb 100644 --- a/lidar_pipeline/rendering.py +++ b/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 mpl_toolkits.axes_grid1.inset_locator import inset_axes +from matplotlib.patches import Polygon as MplPolygon, Rectangle as RectPatch, FancyBboxPatch 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': { - 'cmap': 'coolwarm', - 'title': 'Anomalies Statistiques (MSRM multi-échelle + Moran\'s I)', - '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)', - 'description': 'Détecte uniquement les anomalies statistiquement significatives — filtre le bruit de fond', - 'vmin_mode': 'symmetric', 'sym_pct': (5, 95), + 'svf': { + 'cmap': 'gray_r', + 'title': 'Sky-View Factor (fraction de ciel visible)', + '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étection de micro-relief — fossés sombres, levées claires, complémentaire de l\'openness', + '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_ax = inset_axes(ax, width="5%", height="9%", loc='upper right', - bbox_to_anchor=(-0.03, 0.08, 1, 1), bbox_transform=ax.transAxes) + # North arrow — compass rose style, positioned above the colorbar + compass_bottom = data_bottom + data_height_frac + 0.02 + 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' diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 63e1040..07991dc 100644 --- a/lidar_pipeline/visualizations.py +++ b/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) @@ -1275,4 +1277,180 @@ def generate_flow(dem_file, basename, vis_dir, resolution, shared=None): return output 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 \ No newline at end of file