"""Rendering module: colormap registry, GeoTIFF-to-image conversion, and PDF report generation. Contains: - COLORMAPS: registry mapping filename keywords to (cmap, title, legend, description) - tif_to_png(): convert a GeoTIFF to a WebP/AVIF visualization with legend, scale bar, north arrow - generate_pdf_report(): generate an A3 PDF report with all visualizations """ import logging import time from datetime import datetime from pathlib import Path import numpy as np import rasterio from PIL import Image as PILImage try: from rasterio.warp import transform as warp_transform HAS_WARP = True 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, FancyBboxPatch from matplotlib.colors import ListedColormap try: from cmcrameri import cm as cmc HAS_CMCRAmeri = True except ImportError: HAS_CMCRAmeri = False rcParams['figure.dpi'] = 150 rcParams['savefig.dpi'] = 300 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 # ============================================================ # Each entry: keyword → (cmap, title, legend_label, description, vmin_mode, vmax_mode) # vmin_mode/vmax_mode: 'percentile_X_Y' or '0_max_X' or 'symmetric_X_Y' or 'fixed_0_1' # For RGB images (ortho/topo), special handling is done in tif_to_png. COLORMAPS = { # === Famille RELIEF : rouge=surélévation, bleu=dépression === # Roma (Crameri): perceptually uniform, CVD-friendly, dark center → near-zero values visible # Falls back to RdBu_r if cmcrameri unavailable 'curvature': { 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', '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': 'roma' if HAS_CMCRAmeri else 'RdBu_r', 'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)', 'legend': 'Relief combiné multi-échelles\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nDétecte du micro au macro', 'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, 'lrm': { 'cmap': 'roma' if HAS_CMCRAmeri else 'RdBu_r', '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': 'roma' if HAS_CMCRAmeri else 'RdBu_r', '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': 'roma' if HAS_CMCRAmeri else 'RdBu_r', 'title': 'SAILORE - LRM Auto-Adaptatif', 'legend': 'Relief local adaptatif\nRouge = Surélévation | Bleu = Dépression\n\nLRM = noyau fixe 15m\nMSRM = 5 noyaux fixes\nSAILORE = noyau adapté à la pente\nPlat=grand noyau | Pente=petit noyau', 'description': 'Noyau qui s\'adapte à la pente locale — terrain plat=grand noyau, pente=petit noyau', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, 'aniso_open': { 'cmap': 'roma' if HAS_CMCRAmeri else '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), }, # === Famille OUVERTURE : clair=ouvert, sombre=fermé === 'positive_openness': { 'cmap': 'YlOrBr', 'title': 'Openness Positive (ouverture vers le haut)', 'legend': 'Angle d\'ouverture vers le haut (deg)\nClair = Vue dégagée vers le ciel (sommets, plateaux)\nSombre = Vue bloquée (vallées encaissées)', 'description': 'Ray-tracing 8 directions — complémentaire de la négative pour détecter crêtes', 'vmin_mode': 'percentile', 'vmin_pct': 10, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, 'negative_openness': { 'cmap': 'PuBu', 'title': 'Openness Negative (ouverture vers le bas)', 'legend': 'Angle d\'ouverture vers le bas (deg)\nClair = Surplomb (bords de fossé, grottes)\nSombre = Terrain plat (fonds de vallée)\nMeilleur détecteur de cavités et dolines', 'description': 'Ray-tracing 8 directions — détecte fossés, dolines, souterrains', 'vmin_mode': 'percentile', 'vmin_pct': 10, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, 'svf': { 'cmap': 'bone_r', 'title': 'Sky-View Factor (fraction de ciel visible)', 'legend': 'Proportion de ciel visible depuis chaque point\nClair = Ciel dégagé (sommet, plateau, levée)\nSombre = Ciel masqué (vallée, fossé, tranchée)\nContraste adapté aux valeurs réelles (percentiles 2-98)', 'description': 'Détection de micro-relief — fossés sombres, levées claires, complémentaire de l\'openness', 'vmin_mode': 'percentile', 'vmin_pct': 2, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, # === Famille SCALAIRE : propriétés non divergentes === 'hillshade': { 'cmap': 'gray', 'title': 'Hillshade Multidirectionnel', 'legend': 'Illumination combinée de 6 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre', 'description': 'Ombres portées révélant micro-relief (murs, fossés, terrasses)', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, 'slope': { 'cmap': 'inferno', 'title': 'Pente (Inclinaison du terrain)', 'legend': 'Inclinaison en degrés\nMin: {vmin:.1f}° | Max: {vmax:.1f}°\nClair = Forte pente | Sombre = Terrain plat', 'description': 'Murs, talus et bords ressortent en clair — terrain plat en sombre', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'percentile', 'vmax_pct': 95, }, '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': 'magma', 'title': 'Rugosité Multi-Échelle (3m + 15m)', 'legend': 'Irrégularité du terrain combinée fine + large\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nCombine rugosité fine 3m (70%) + large 15m (30%)', 'description': 'Mesure la variabilité locale — surfaces anthropiques lisses vs naturelles rugueuses', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'percentile', 'vmax_pct': 97, }, 'wavelet': { 'cmap': 'cividis', 'title': 'Ondelette Mexican Hat (CWT 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', 'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, } # RGB entries (ortho/topo) are handled specially RGB_LEGENDS = { 'ortho': { 'title': 'Photographie Aérienne IGN', 'legend': 'Orthophotographie\nImage aérienne', 'description': 'Photographie aérienne IGN (Orthophoto)', }, 'topo': { 'title': 'Carte Topographique IGN', 'legend': 'Carte IGN\nPlan topographique', 'description': 'Carte topographique IGN (Plan IGN)', }, } def _apply_colormap(data, tif_file): """Apply the registered colormap normalization to data based on filename. Returns (data, cmap, title, legend_label, description, is_rgb). """ name = str(tif_file).lower() # Check for RGB first for key in RGB_LEGENDS: if key in name: info = RGB_LEGENDS[key] return data, None, info['title'], info['legend'], info['description'], True # Find matching colormap — sort by key length descending so 'mslrm' matches before 'lrm' for key in sorted(COLORMAPS.keys(), key=len, reverse=True): info = COLORMAPS[key] if key in name: valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten()) valid_data = valid_data[~np.isnan(valid_data)] if len(valid_data) == 0: logger.warning(f" Aucune donnée valide pour {Path(tif_file).name} — colormap ignorée") return data, 'terrain', Path(tif_file).stem.replace('_', ' ').title(), '', '', False vmin = vmax = None # Compute vmin/vmax based on mode if info['vmin_mode'] == 'fixed': vmin = info['vmin_val'] elif info['vmin_mode'] == 'percentile': vmin = np.percentile(valid_data, info['vmin_pct']) elif info['vmin_mode'] == 'symmetric': vmax_abs = max(abs(np.percentile(valid_data, info['sym_pct'][0])), abs(np.percentile(valid_data, info['sym_pct'][1])), 0.001) vmin = -vmax_abs vmax = vmax_abs # symmetric mode sets both vmin and vmax if vmax is None: # Only compute vmax if not already set by symmetric mode if info.get('vmax_mode') == 'fixed': vmax = info['vmax_val'] elif info.get('vmax_mode') == 'percentile': vmax = np.percentile(valid_data, info['vmax_pct']) elif info.get('vmax_mode') == 'symmetric': vmax_abs = max(abs(np.percentile(valid_data, info['sym_pct'][0])), abs(np.percentile(valid_data, info['sym_pct'][1])), 0.001) vmax = vmax_abs # Apply normalization if vmin is not None and vmax is not None: data = np.clip((data - vmin) / max(vmax - vmin, 0.001), 0, 1) legend = info['legend'].format(vmin=vmin or 0, vmax=vmax or 0) return data, info['cmap'], info['title'], legend, info['description'], False # Default: terrain colormap with percentile stretch valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten()) valid_data = valid_data[~np.isnan(valid_data)] if len(valid_data) == 0: return data, 'terrain', Path(tif_file).stem.replace('_', ' ').title(), '', '', False p2, p98 = np.percentile(valid_data, (2, 98)) data = np.clip((data - p2) / (p98 - p2), 0, 1) title = Path(tif_file).stem.replace('_', ' ').title() 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 ~30km around the processed zone, giving regional 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: Tuple (image_array, bounds_dict) where bounds_dict has keys 'min_x', 'max_x', 'min_y', 'max_y' in Lambert 93, 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 zoom 10 for context (~150m/px — wide view, fast download) context_zoom = 10 # Expand bounds to ~30km for regional context # Large enough to see surrounding towns/rivers, small enough that # the processed zone (typically 1km) is clearly visible as a red rectangle context_half = 15000 # 15km each side = 30km total context_min_x = center_x - context_half context_max_x = center_x + context_half context_min_y = center_y - context_half context_max_y = center_y + context_half result = download_ign_tiles( context_min_x, context_max_x, context_min_y, context_max_y, layer='GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2', zoom_level=context_zoom, min_zoom=8 ) if result is not None: bounds = { 'min_x': context_min_x, 'max_x': context_max_x, 'min_y': context_min_y, 'max_y': context_max_y, } cached = (result, bounds) _location_map_cache[cache_key] = cached return cached 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. Returns (scale_m, label) where label is like '100 m' or '500 m' or '1 km'. """ nice_scales = [50, 100, 200, 500, 1000, 2000, 5000, 10000] # Pick the largest scale <= 20% of extent chosen = nice_scales[0] for s in nice_scales: if s <= extent_m * 0.20: chosen = s if chosen >= 1000: return chosen, f"{chosen // 1000} km" return chosen, f"{chosen} m" def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, quality=98, output_format='avif'): """Convert GeoTIFF to visualization image (WebP or AVIF) with GPS coordinates, legend, and scale bar. Args: tif_file: Path to input GeoTIFF. vis_dir: Output directory for the image file. resolution: Grid resolution in m/px. keep_tif: If True, keep the source TIFF after conversion. source_info: Dict with method/date/basename for metadata. quality: Image quality (1-100). Use 100 for lossless. Default 95. output_format: Output format ('webp' or 'avif'). Default 'webp'. Returns: Path to output image file, or None on failure. """ if not tif_file or not tif_file.exists(): return None ext = 'avif' if output_format == 'avif' else 'webp' output_file = vis_dir / f"{tif_file.stem}.{ext}" try: with rasterio.open(tif_file) as src: is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo')) if is_rgb: data = src.read([1, 2, 3]) data = np.moveaxis(data, 0, -1) else: data = src.read(1) nodata = src.nodata transform = src.transform crs = src.crs if is_rgb: height, width, _ = data.shape else: height, width = data.shape top_left_x = transform.c top_left_y = transform.f pixel_size_x = transform.a pixel_size_y = abs(transform.e) min_x = top_left_x max_x = top_left_x + width * pixel_size_x max_y = top_left_y min_y = top_left_y - height * pixel_size_y # GPS coordinates gps_coords = {} if HAS_WARP and crs is not None: try: l93_xs = [min_x, max_x, min_x, max_x] l93_ys = [max_y, max_y, min_y, min_y] lons, lats = warp_transform(crs, 'EPSG:4326', l93_xs, l93_ys) gps_coords = { 'NW': (lats[0], lons[0]), 'NE': (lats[1], lons[1]), 'SW': (lats[2], lons[2]), 'SE': (lats[3], lons[3]), } n_ticks = 5 tick_l93_x = np.linspace(min_x, max_x, n_ticks) tick_l93_y_bottom = np.full(n_ticks, min_y) tick_lons, tick_lats = warp_transform(crs, 'EPSG:4326', tick_l93_x, tick_l93_y_bottom) gps_coords['x_ticks'] = list(zip(tick_lons, tick_lats)) tick_l93_y = np.linspace(min_y, max_y, n_ticks) tick_l93_x_left = np.full(n_ticks, min_x) tick_lons_y, tick_lats_y = warp_transform(crs, 'EPSG:4326', tick_l93_x_left, tick_l93_y) gps_coords['y_ticks'] = list(zip(tick_lons_y, tick_lats_y)) except Exception: gps_coords = {} if nodata is not None and not is_rgb: data = np.ma.masked_where((data == nodata) | np.isnan(data), data) if not is_rgb: valid_data = np.asarray(data.compressed() if hasattr(data, 'compressed') else data.flatten()) valid_data = valid_data[~np.isnan(valid_data)] # Track NaN mask before converting to plain ndarray nan_mask = None if not is_rgb: if isinstance(data, np.ma.MaskedArray): nan_mask = data.mask.copy() data = np.ma.filled(data, np.nan) elif np.any(np.isnan(data)): nan_mask = np.isnan(data) # For rendering: replace NaN with neutral value to avoid interpolation halos if nan_mask is not None and np.any(nan_mask) and len(valid_data) > 0: fill_value = float(np.median(valid_data)) data[nan_mask] = fill_value nan_mask = nan_mask # keep for later # Apply colormap data, cmap, title, legend_label, description, is_rgb_result = _apply_colormap(data, tif_file) # Apply NaN mask: make zones without data transparent has_nan_mask = nan_mask is not None and not is_rgb_result if has_nan_mask: # data is normalized 0-1 from _apply_colormap; apply cmap to get RGBA # Save the colormap for colorbar before converting to RGBA saved_cmap = plt.get_cmap(cmap) if isinstance(cmap, str) else cmap saved_vmin = float(np.nanmin(data)) if not nan_mask.all() else 0 saved_vmax = float(np.nanmax(data)) if not nan_mask.all() else 1 rgba = saved_cmap(data) # (H, W, 4) float RGBA rgba[nan_mask, 3] = 0.0 # transparent where no data data = rgba is_rgba = True else: is_rgba = False saved_cmap = None saved_vmin = None saved_vmax = None # Create figure with FIXED layout for consistent data area position # All visualizations use the same axes positions so they can be overlaid fig_width = max(20, width / 150) fig_width = min(fig_width, 40) fig_height = fig_width * 0.7 + 2.0 # Fixed header + footer space fig = plt.figure(figsize=(fig_width, fig_height), facecolor='white') # Fixed data area position — identical for ALL visualization types # This ensures overlay/superposition works across all WebP images data_left = 0.08 data_bottom = 0.19 data_width_frac = 0.74 data_height_frac = 0.71 ax = fig.add_axes([data_left, data_bottom, data_width_frac, data_height_frac]) if is_rgba or is_rgb: im = ax.imshow(data, aspect='equal', origin='upper', interpolation='bilinear') else: im = ax.imshow(data, cmap=cmap, aspect='equal', origin='upper', interpolation='bilinear') ax.set_title(f"{title}\n{description}", fontsize=15, fontweight='bold', pad=10) # 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, 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, fontsize=9, fontweight='bold', rotation=90, verticalalignment='center', horizontalalignment='center', wrap=True) cbar_ax.set_frame_on(False) elif is_rgba and saved_cmap is not None: 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([]) cbar = plt.colorbar(sm, cax=cbar_ax) cbar.ax.tick_params(labelsize=9, width=1.5) cbar.outline.set_linewidth(1.5) cbar.set_label(legend_label, fontsize=10, fontweight='bold') else: 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) cbar.set_label(legend_label, fontsize=10, fontweight='bold') # GPS coordinate ticks if gps_coords and 'x_ticks' in gps_coords: x_pixel_positions = np.linspace(0, width - 1, len(gps_coords['x_ticks'])) x_labels = [f"{lon:.5f}E" for lon, lat in gps_coords['x_ticks']] ax.set_xticks(x_pixel_positions) ax.set_xticklabels(x_labels, fontsize=7, rotation=30) ax.set_xlabel('Longitude', fontsize=9, fontweight='bold') y_pixel_positions = np.linspace(0, height - 1, len(gps_coords['y_ticks'])) y_labels = [f"{lat:.5f}N" for lon, lat in gps_coords['y_ticks']] ax.set_yticks(y_pixel_positions) ax.set_yticklabels(y_labels, fontsize=7) ax.set_ylabel('Latitude', fontsize=9, fontweight='bold') else: x_ticks_count = 5 x_positions = np.linspace(0, width - 1, x_ticks_count) x_labels = [f"{(min_x + xp * pixel_size_x)/1000:.1f}" for xp in x_positions] ax.set_xticks(x_positions) ax.set_xticklabels(x_labels, fontsize=8) ax.set_xlabel('Est (km) - Lambert 93', fontsize=9, fontweight='bold') y_ticks_count = 5 y_positions = np.linspace(0, height - 1, y_ticks_count) y_labels = [f"{(max_y - yp * pixel_size_y)/1000:.1f}" for yp in y_positions] ax.set_yticks(y_positions) ax.set_yticklabels(y_labels, fontsize=8) ax.set_ylabel('Nord (km) - Lambert 93', fontsize=9, fontweight='bold') ax.tick_params(axis='both', which='both', direction='out', length=3, width=0.8, colors='black') for spine in ax.spines.values(): spine.set_visible(True) spine.set_color('black') spine.set_linewidth(0.8) # 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') north_ax.plot([0, 0], [0.0, 1.0], color='#b22222', linewidth=2.0, zorder=10) north_ax.add_patch(MplPolygon([[0, 0.3], [-0.2, 0.7], [0, 1.0], [0.2, 0.7]], closed=True, facecolor='#b22222', edgecolor='#b22222', zorder=9)) # Cardinal ticks for angle, label in [(90, ''), (0, 'E'), (180, 'O'), (270, 'S')]: rad = np.radians(angle) r_text = 1.25 north_ax.plot([0.85*np.cos(rad), 1.05*np.cos(rad)], [0.85*np.sin(rad), 1.05*np.sin(rad)], color='#555555', linewidth=0.8, zorder=5) if label: north_ax.text(r_text*np.cos(rad), r_text*np.sin(rad), label, fontsize=7, ha='center', va='center', color='#555555') # Bottom info bar — enriched with source, method, date info_ax = fig.add_axes([data_left, 0.015, data_width_frac + cbar_width + 0.02, 0.09]) info_ax.axis('off') extent_km_x = (max_x - min_x) / 1000 extent_km_y = (max_y - min_y) / 1000 if is_rgb: alt_min = alt_max = 0 else: alt_min = float(np.nanmin(valid_data)) if len(valid_data) > 0 else 0 alt_max = float(np.nanmax(valid_data)) if len(valid_data) > 0 else 0 # Build info lines line1_parts = [] if gps_coords: nw_lat, nw_lon = gps_coords['NW'] se_lat, se_lon = gps_coords['SE'] line1_parts.append(f"GPS: {nw_lat:.5f}°N {nw_lon:.5f}°E — {se_lat:.5f}°N {se_lon:.5f}°E") else: line1_parts.append(f"X: {min_x:.0f}–{max_x:.0f} Y: {min_y:.0f}–{max_y:.0f}") line1_parts.append(f"EPSG:2154") line1_parts.append(f"Res: {resolution}m/px") line1_parts.append(f"Emprise: {extent_km_x:.1f}×{extent_km_y:.1f}km") if not is_rgb: line1_parts.append(f"Alt: {alt_min:.1f}–{alt_max:.1f}m") line2_parts = [] line2_parts.append("Source: LiDAR HD IGN") if source_info: if source_info.get('method'): line2_parts.append(f"Classif.: {source_info['method'].upper()}") if source_info.get('date'): line2_parts.append(f"Date: {source_info['date']}") else: line2_parts.append(datetime.now().strftime("Date: %Y-%m-%d")) info_text_line1 = " | ".join(line1_parts) info_text_line2 = " | ".join(line2_parts) info_ax.text(0.01, 0.7, info_text_line1, transform=info_ax.transAxes, fontsize=8, verticalalignment='center', family='monospace', bbox=dict(boxstyle='round,pad=0.2', facecolor='#f0f0f0', edgecolor='#aaaaaa', alpha=0.95)) info_ax.text(0.01, 0.2, info_text_line2, transform=info_ax.transAxes, fontsize=7.5, verticalalignment='center', family='monospace', color='#444444', bbox=dict(boxstyle='round,pad=0.2', facecolor='#f8f8f8', edgecolor='#cccccc', alpha=0.9)) # Scale bar — adaptive with alternating black/white segments # Position: left of the location map to avoid overlap extent_m_x = max_x - min_x scale_m, scale_label = _nice_scale(extent_m_x) pixels_per_meter = 1.0 / pixel_size_x scale_px = int(scale_m * pixels_per_meter) n_segments = 5 segment_px = scale_px / n_segments bar_bottom_y = 0.55 bar_top_y = 0.85 bar_height = bar_top_y - bar_bottom_y # Place scale bar so it ends before the location map (map starts at x=0.82 in fig coords) # map_ax occupies [0.82, 0.02, 0.16, 0.13] in figure coords # info_ax occupies [data_left, 0.015, width, 0.09] # Scale bar end in fig coords = info_ax.left + (scale_start_x + scale_px/width) * info_ax.width # We need: info_ax.left + (scale_start_x + scale_px/width) * info_ax.width < 0.80 scale_end_frac = scale_px / width # fraction of info_ax width info_ax_width = data_width_frac + cbar_width + 0.02 # Calculate scale_start_x so scale bar ends at fig_x = 0.78 (leaving gap before map at 0.82) max_scale_end_fig = 0.78 scale_end_in_info = (max_scale_end_fig - data_left) / info_ax_width scale_start_x = max(0.05, scale_end_in_info - scale_end_frac) for seg_i in range(n_segments): color = 'black' if seg_i % 2 == 0 else 'white' seg_left = scale_start_x + seg_i * segment_px / width seg_width_frac = segment_px / width info_ax.add_patch(RectPatch((seg_left, bar_bottom_y), seg_width_frac, bar_height, facecolor=color, edgecolor='black', linewidth=0.5, transform=info_ax.transAxes, clip_on=False)) info_ax.text(scale_start_x + scale_px / (2 * width), bar_top_y + 0.12, f"{scale_label}", ha='center', va='bottom', fontsize=8, fontweight='bold', transform=info_ax.transAxes) # Scale end ticks info_ax.plot([scale_start_x, scale_start_x], [bar_bottom_y - 0.05, bar_top_y + 0.05], color='black', linewidth=1, transform=info_ax.transAxes, clip_on=False) info_ax.plot([scale_start_x + scale_px / width, scale_start_x + scale_px / width], [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 # Positioned in lower-right corner, above the info bar map_ax = fig.add_axes([0.82, 0.02, 0.16, 0.13]) # Try to download a wide-area IGN topo map for location context location_result = _download_location_map(min_x, max_x, min_y, max_y) if location_result is not None: location_map, loc_bounds = location_result # Draw IGN topo map as background with correct bounds map_ax.imshow(location_map, aspect='equal', extent=[ loc_bounds['min_x'], loc_bounds['max_x'], loc_bounds['min_y'], loc_bounds['max_y'] ]) # 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' save_dpi = 200 if width > 3000 else 150 png_temp = vis_dir / f"{tif_file.stem}_temp.png" try: plt.savefig(png_temp, dpi=save_dpi, facecolor='white', format='png') finally: plt.close() img = PILImage.open(str(png_temp)) pil_format = 'AVIF' if output_format == 'avif' else 'WEBP' if quality >= 100: img.save(str(output_file), format=pil_format, lossless=True) else: img.save(str(output_file), format=pil_format, quality=quality) png_temp.unlink(missing_ok=True) # Delete source TIFF (unless --keep-tif) if not keep_tif: tif_file.unlink(missing_ok=True) return output_file except Exception as e: logger.error(f" Erreur conversion {ext.upper()}: {e}", exc_info=True) return None def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): """Generate A3 PDF report for a LiDAR file with all visualizations. Page 1: Mise en situation (ortho + topo IGN side by side) Pages 2+: Other visualizations (2 per page) Args: basename: Base name for the report file. vis_dir: Directory containing WebP visualization files. pdf_dir: Directory for output PDF. resolution: Grid resolution (used in info text). Returns: Path to PDF file, or None on failure. """ from matplotlib.backends.backend_pdf import PdfPages pdf_file = pdf_dir / f"{basename}_rapport.pdf" logger.info(f" → Génération rapport PDF A3: {pdf_file.name}") t0 = time.time() # Look for WebPs in per-file subdirectory first, then fallback to main dir file_vis_dir = vis_dir / basename if file_vis_dir.exists(): png_files = sorted(file_vis_dir.glob("*.webp")) else: png_files = sorted(vis_dir.glob(f"{basename}_*.webp")) if not png_files: logger.warning(f" ✗ Aucune image trouvée pour {basename}") return None # Categorize situ_files = [] analysis_files = [] for f in png_files: name = f.stem.lower() if 'ortho' in name: situ_files.insert(0, f) elif 'topo' in name: situ_files.append(f) else: analysis_files.append(f) # Sort analysis files by archaeological priority order = ['mslrm', 'svf', 'negative_openness', 'positive_openness', 'aniso_open', 'sailore', 'hillshade_multi', 'lrm', 'tpi', 'slope', 'curvature', 'aspect', 'roughness', 'wavelet'] def sort_key(f): name = f.stem.lower() for i, key in enumerate(order): if key in name: return i return len(order) analysis_files.sort(key=sort_key) a3_w, a3_h = 16.54, 11.69 try: with PdfPages(str(pdf_file)) as pdf: # Page 1: Mise en situation if situ_files: fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white') n_situ = len(situ_files) if n_situ == 2: gs = fig.add_gridspec(1, 2, wspace=0.05, left=0.03, right=0.97, top=0.92, bottom=0.06) else: gs = fig.add_gridspec(1, max(n_situ, 1), wspace=0.05, left=0.03, right=0.97, top=0.92, bottom=0.06) fig.text(0.5, 0.97, f"Mise en situation - {basename}", fontsize=20, fontweight='bold', ha='center', va='top') for i, f in enumerate(situ_files): ax = fig.add_subplot(gs[0, i]) img = np.array(PILImage.open(str(f)).convert('RGB')) ax.imshow(img) ax.axis('off') title = f.stem.replace(basename + '_', '').replace('_', ' ').title() ax.set_title(title, fontsize=12, fontweight='bold', pad=5) pdf.savefig(fig, dpi=150) plt.close(fig) # Pages 2+: Analysis maps (2 per page) for page_start in range(0, len(analysis_files), 2): page_files = analysis_files[page_start:page_start + 2] fig = plt.figure(figsize=(a3_w, a3_h), facecolor='white') if len(page_files) == 2: gs = fig.add_gridspec(1, 2, wspace=0.08, left=0.03, right=0.97, top=0.93, bottom=0.05) else: gs = fig.add_gridspec(1, 1, left=0.05, right=0.95, top=0.93, bottom=0.05) for i, f in enumerate(page_files): ax = fig.add_subplot(gs[0, i]) img = np.array(PILImage.open(str(f)).convert('RGB')) ax.imshow(img) ax.axis('off') title = f.stem.replace(basename + '_', '').replace('_', ' ').title() ax.set_title(title, fontsize=11, fontweight='bold', pad=3) page_num = (page_start // 2) + 2 fig.text(0.99, 0.01, f"Page {page_num}", fontsize=8, ha='right', va='bottom', color='gray') pdf.savefig(fig, dpi=150) plt.close(fig) logger.info(f" ✓ Rapport PDF terminé ({time.time()-t0:.1f}s)") return pdf_file except Exception as e: logger.error(f" ✗ Erreur PDF: {e}", exc_info=True) return None