"""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 from matplotlib.ticker import ScalarFormatter 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, 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. # Les textes (title/legend/description) viennent de VIZ_LEGENDS (index.py, # source unique partagée avec l'export multi-dalles) — fusion ci-dessous. COLORMAPS = { # === Famille RELIEF : rouge=surélévation, bleu=dépression === # Diverging: rouge vif=positif, bleu vif=négatif, blanc=plat 'mslrm': { 'cmap': 'seismic', 'vmin_mode': 'fixed', 'vmin_val': -3, 'vmax_mode': 'fixed', 'vmax_val': 3, }, 'sailore': { 'cmap': 'seismic', 'vmin_mode': 'fixed', 'vmin_val': -3, 'vmax_mode': 'fixed', 'vmax_val': 3, }, # === Famille OUVERTURE : séquentiel, valeurs normalisées === 'positive_openness': { 'cmap': 'YlOrBr', 'vmin_mode': 'fixed', 'vmin_val': -3, 'vmax_mode': 'fixed', 'vmax_val': 3, }, 'negative_openness': { 'cmap': 'PuBu', 'vmin_mode': 'fixed', 'vmin_val': -3, 'vmax_mode': 'fixed', 'vmax_val': 3, }, 'svf': { 'cmap': 'hot_r', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, # === Famille SCALAIRE : propriétés non divergentes === # Clé = mot-clé exact du nom de fichier de sortie (hillshade_multi). 'hillshade_multi': { 'cmap': 'gray', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, 'slope': { 'cmap': 'inferno', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 30, }, 'aspect': { 'cmap': 'twilight', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 360, }, 'roughness': { 'cmap': 'plasma', # vmax fixe (p98 médian mesuré sur 20 dalles réelles) : un vmax au # percentile par tuile rendait l'échelle non jointive entre tuiles. 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 3.8, }, 'wavelet': { 'cmap': 'inferno', # Nœuds (percentile → valeur) mesurés sur 20 tuiles réelles par # résolution avec l'algorithme actuel (détendage 35 m, échelles # 1-50 m) : médiane inter-tuiles des percentiles par tuile, robuste # aux tuiles très structurées. Distribution plus large qu'avant # détendage : le fond macro-relief supprimé, les petites structures # percent bien plus au-dessus du bruit (p98 ≈ 9 vs 1,65 avant). 'knots': { 0.5: ([0.089, 0.213, 0.296, 0.445, 0.602, 0.783, 1.0, 1.274, 1.661, 2.32, 3.842, 5.661, 8.936, 11.93, 15.06], [0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.98, 0.99, 0.995]), 0.2: ([0.088, 0.213, 0.296, 0.445, 0.602, 0.783, 1.0, 1.287, 1.694, 2.358, 3.846, 5.687, 8.971, 11.928, 14.971], [0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.98, 0.99, 0.995]), }, }, 'flow_acc': { 'cmap': 'YlGn', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, 'anomaly': { 'cmap': 'YlOrRd', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, 'solar': { 'cmap': 'gray', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, } # RGB entries (ortho/topo) are handled specially RGB_LEGENDS = { 'ortho': {}, 'topo': {}, } # Fusion des textes de légende (titre / lecture du rendu / méthode de calcul) # depuis la source unique VIZ_LEGENDS (index.py, sans dépendance lourde). from .index import VIZ_LEGENDS for _key, _info in COLORMAPS.items(): _info.update(VIZ_LEGENDS[_key]) for _key, _info in RGB_LEGENDS.items(): _info.update(VIZ_LEGENDS[_key]) def _apply_colormap(data, tif_file, resolution=None): """Apply the registered colormap normalization to data based on filename. Returns (data, cmap, title, legend_label, description, is_rgb, vmin, vmax) où vmin/vmax sont les bornes physiques de la plage rendue (None si sans objet). """ 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, None, None # Find matching colormap — tri par longueur décroissante : en cas de # chevauchement de mots-clés dans un nom de fichier, le plus long prime 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, None, None vmin = vmax = None knots = info.get('knots') if knots is not None: # Étalonnage quantile figé (appariement d'histogramme global, # cf. normalisation radiométrique des mosaïques) : fonction de # transfert en nœuds mesurés une fois sur un échantillon de # tuiles — même valeur → même couleur sur toutes les tuiles, # toute la palette utilisée, insensible aux queues locales. if isinstance(knots, dict): key = min(knots, key=lambda r: abs(float(r) - float(resolution or 0.5))) knots = knots[key] kv, kt = knots data = np.interp(np.asarray(data, dtype=float), kv, kt, left=0.0, right=1.0) vmin, vmax = kv[0], kv[-1] else: # 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, vmin, vmax # 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)) # Garde contre les tuiles quasi constantes : plage nulle → division par zéro span = max(p98 - p2, 1e-6) data = np.clip((data - p2) / span, 0, 1) title = Path(tif_file).stem.replace('_', ' ').title() return data, 'terrain', title, 'Altitude normalisée', '', False, p2, p98 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. """ # 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) # Éviction FIFO : chaque entrée ~12 Mo (carte ~30 km au zoom 10) — # un run long sur plusieurs zones ne doit pas empiler sans limite while len(_location_map_cache) >= 4: _location_map_cache.pop(next(iter(_location_map_cache))) _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 = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000] # Pick the largest scale that fits within 20% of extent max_scale = extent_m * 0.20 chosen = nice_scales[0] for s in nice_scales: if s <= max_scale: chosen = s else: break 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 98. output_format: Output format ('webp' or 'avif'). Default 'avif'. 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, cmap_vmin, cmap_vmax = _apply_colormap(data, tif_file, resolution=resolution) # 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 # Bornes physiques de la colorbar (unités réelles, pas 0-1) if cmap_vmin is not None and cmap_vmax is not None: saved_vmin, saved_vmax = float(cmap_vmin), float(cmap_vmax) else: 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}", fontsize=14, fontweight='bold', pad=10) if description: ax.text(0.5, 1.04, description, transform=ax.transAxes, fontsize=10, fontstyle='italic', color='#555555', ha='center', va='bottom') # Colorbar/legend area — full height alongside data cbar_left = data_left + data_width_frac + 0.02 cbar_width = 0.04 cbar_bottom = data_bottom cbar_height = data_height_frac 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.ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False)) 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.ax.yaxis.set_major_formatter(ScalarFormatter(useOffset=False)) 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 in bottom-right corner of data area # Semi-transparent background for readability over any data north_ax = fig.add_axes([data_left + data_width_frac - 0.07, data_bottom + 0.01, 0.06, 0.14], facecolor='none') north_ax.set_xlim(-1.5, 1.5) north_ax.set_ylim(-1.5, 1.5) north_ax.axis('off') north_ax.set_aspect('equal') # Compass rose centered at (0, 0) — all 4 cardinals equidistant from center # Semi-transparent white background circle circle_bg = plt.Circle((0, 0), 1.0, facecolor='white', edgecolor='#888888', linewidth=0.5, alpha=0.7, zorder=1) north_ax.add_patch(circle_bg) # N arrow (pointing up = North) north_ax.annotate('N', xy=(0, 1.35), fontsize=9, fontweight='bold', ha='center', va='bottom', color='#b22222', zorder=10) north_ax.plot([0, 0], [-0.5, 1.0], color='#b22222', linewidth=2.0, zorder=10) north_ax.add_patch(MplPolygon([[0, 0.5], [-0.2, 0.7], [0, 1.0], [0.2, 0.7]], closed=True, facecolor='#b22222', edgecolor='#b22222', zorder=9)) # Cardinal ticks — all centered at (0, 0) for angle, label in [(90, 'N'), (0, 'E'), (180, 'O'), (270, 'S')]: rad = np.radians(angle) north_ax.plot([1.0*np.cos(rad), 1.2*np.cos(rad)], [1.0*np.sin(rad), 1.2*np.sin(rad)], color='#555555', linewidth=0.8, zorder=5) if label: north_ax.text(1.35*np.cos(rad), 1.35*np.sin(rad), label, fontsize=6, ha='center', va='center', color='#555555', zorder=5) # 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") # Round resolution to avoid ugly decimals like 0.499999959 res_display = round(resolution, 2) if resolution < 1 else round(resolution, 1) line1_parts.append(f"Res: {res_display}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 figure to in-memory buffer (avoids disk I/O of temp PNG) save_dpi = 200 if width > 3000 else 150 from io import BytesIO buf = BytesIO() try: plt.savefig(buf, dpi=save_dpi, facecolor='white', format='png') finally: plt.close() buf.seek(0) img = PILImage.open(buf) 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) # 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 tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=98, output_format='avif'): """Convert GeoTIFF to a cropped visualization image (no legend, no overlay). Applies colormap and saves the image as a pure 1×1 km square. Used for grid/map display where images must tile seamlessly. 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. quality: Image quality (1-100). Use 100 for lossless. output_format: Output format ('webp' or 'avif'). 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 → NaN : sinon les pixels bord (ex. -9999, 3.4e38) # polluent les percentiles d'étalonnage de la colormap if src.nodata is not None: data = data.astype(np.float32, copy=False) data[data == src.nodata] = np.nan # Apply colormap normalization data, cmap_name, title, legend_label, description, is_rgb_result, _cvmin, _cvmax = _apply_colormap(data, tif_file, resolution=resolution) if not is_rgb_result: data = np.where(np.isnan(data), 0.0, data) # Convert to RGB using colormap if is_rgb_result: # RGB images are already in RGB (uint8 depuis le TIF IGN, ou float 0-1) if data.dtype == np.uint8: rgb_data = data else: rgb_data = (np.clip(data, 0, 1) * 255).astype(np.uint8) else: # Normalize data to 0-1 range for colormap cmap = plt.get_cmap(cmap_name) rgb_float = cmap(data.clip(0, 1)) rgb_data = (rgb_float[:, :, :3] * 255).astype(np.uint8) # Save as AVIF/WebP img = PILImage.fromarray(rgb_data) 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) # 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 crop {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 images in per-file subdirectory first, then fallback to main dir file_vis_dir = vis_dir / basename png_files = [] if file_vis_dir.exists(): png_files = sorted(file_vis_dir.glob("*.avif")) + sorted(file_vis_dir.glob("*.webp")) else: png_files = sorted(vis_dir.glob(f"{basename}_*.avif")) + sorted(vis_dir.glob(f"{basename}_*.webp")) # Deduplicate in case both formats exist seen = set() unique_files = [] for f in png_files: if f not in seen: seen.add(f) unique_files.append(f) if not unique_files: logger.warning(f" ✗ Aucune image trouvée pour {basename}") return None png_files = unique_files # 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', 'sailore', 'hillshade_multi', 'flow_acc', 'solar', 'slope', 'roughness', 'wavelet', 'aspect', 'anomaly'] 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]) with PILImage.open(str(f)) as _pf: img = np.array(_pf.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]) with PILImage.open(str(f)) as _pf: img = np.array(_pf.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