"""Rendering module: colormap registry and GeoTIFF-to-image conversion. Contains: - COLORMAPS: registry mapping filename keywords to a normalization (cmap, vmin/vmax or quantile knots) plus the texts merged from VIZ_LEGENDS - tif_to_crop(): convert a GeoTIFF to a bare 1 km AVIF/WebP tile (used by the pipeline for the map mosaic) - tif_to_png(): convert a GeoTIFF to an annotated AVIF/WebP sheet (legend, scale bar, north arrow) — not called by the pipeline - generate_pdf_report(): legacy per-tile A3 PDF report, no longer called (the map's PDF export lives in export_pdf.py) """ 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/relief_oriente), see RGB_LEGENDS below. # The texts (title/legend/description) come from VIZ_LEGENDS (index.py, the # single source shared with the map UI, the TileJSON and the PDF export) — # merged below. COLORMAPS = { # === RELIEF family: red = raised, blue = depression === # Diverging: bright red = positive, bright blue = negative, white = flat '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, }, # === OPENNESS family: sequential, normalized values === '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, }, # === SCALAR family: non-diverging properties === # Key = exact keyword of the output file name (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', # Fixed vmax (median p98 measured on 20 real tiles): a per-tile # percentile vmax made the scale inconsistent from tile to tile. 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 3.8, }, 'wavelet': { 'cmap': 'inferno', # Knots (value → percentile) measured on 20 real tiles per # resolution with the current algorithm (35 m detrending, 1-50 m # scales): cross-tile median of the per-tile percentiles, robust to # highly structured tiles. Wider distribution than before # detrending: with the macro-relief background removed, small # structures stand out far more above the noise (p98 ≈ 9 vs 1.65 before). '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, }, # Ground point density: level 0..15 (visualizations.density_levels) 'densite_sol': { 'cmap': 'gray', 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 15, }, } # Layers made of coded flat gray levels (discrete levels): grayscale image # (mode L) encoded LOSSLESSLY as WebP (tif_to_crop; Pillow ignores # lossless=True in AVIF) — AVIF q55-60 shifts ~20% of the pixels by up to 4 # levels on these flat areas. They stay at their own resolution (1 m for # the density: ~150–260 KB per tile). LOSSLESS_GRAY_KEYWORDS = ('densite_sol',) # RGB entries (ortho/topo/relief_oriente) are handled specially RGB_LEGENDS = { 'ortho': {}, 'topo': {}, 'relief_oriente': {}, # computed RGB (visualizations.generate_relief_oriente) } RGB_KEYWORDS = tuple(RGB_LEGENDS) # AVIF encoding speed (libavif, 0 = slow/compact … 10 = fast). Measured on # a 5000 × 5000 px tile (q60): default 4.1 s; speed 9 0.6 s for +3% size and # −0.3 dB PSNR, invisible. Encoding was the longest step of rendering a # layer. AVIF_SPEED = 9 # Merge the legend texts (title / how to read the rendering / computation # method) from the single source VIZ_LEGENDS (index.py, no heavy dependency). from .index import VIZ_LEGENDS, parse_basename_coords for _key, _info in COLORMAPS.items(): _info.update(VIZ_LEGENDS[_key]) for _key, _info in RGB_LEGENDS.items(): _info.update(VIZ_LEGENDS[_key]) def _core_tile_window(tif_file, src): """Raster window of the nominal 1 km tile (edge stitching). A TIF produced over an extended footprint (edge buffer filled with the neighboring tiles, see --edge-buffer in dtm.py) is cropped to the exact LHD tile: the final images remain 1 km squares aligned on the multi-tile grid, with no artifact from one tile to the next. Returns None if the TIF does not extend beyond it (nothing to crop) or if the name carries no LHD coordinates. """ coords = parse_basename_coords(Path(tif_file).stem) if coords is None: return None col_km, row_km = coords # LHD grid: (col, row) = north-west corner in km → X ∈ [col, col+1] km, # Y ∈ [row-1, row] km (north edge = row). west = float(col_km) * 1000.0 north = float(row_km) * 1000.0 south = north - 1000.0 try: from rasterio.windows import Window, from_bounds win = from_bounds(west, south, west + 1000.0, north, src.transform) win = win.round_offsets().round_lengths().intersection( Window(0, 0, src.width, src.height)) except Exception: return None if win.width < src.width or win.height < src.height: return win return None 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) where vmin/vmax are the physical bounds of the rendered range (None when not applicable). """ 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 — sorted by decreasing length: when keywords # overlap in a file name, the longest one wins 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" No valid data in {Path(tif_file).name} — colormap skipped") return data, 'terrain', Path(tif_file).stem.replace('_', ' ').title(), '', '', False, None, None vmin = vmax = None knots = info.get('knots') if knots is not None: # Frozen quantile calibration (global histogram matching, # as in radiometric normalization of mosaics): transfer # function given as knots measured once on a sample of # tiles — same value → same color on every tile, the whole # palette is used, insensitive to local tails. 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)) # Guard against near-constant tiles: zero range → division by zero 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, 'Normalized elevation', '', 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) # FIFO eviction: each entry ~12 MB (~30 km map at zoom 10) — a long # run over several areas must not pile up without limit 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" IGN location map unavailable: {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=60, 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). Default 60. 100 requests lossless, which only WebP honors (Pillow ignores lossless=True in AVIF). 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 RGB_KEYWORDS) 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 # Edge stitching: crop to the nominal 1 km tile (GPS coordinates, # scale and legend follow the crop). core_win = _core_tile_window(tif_file, src) if core_win is not None: rows = slice(core_win.row_off, core_win.row_off + core_win.height) cols = slice(core_win.col_off, core_win.col_off + core_win.width) data = data[rows, cols, :] if is_rgb else data[rows, cols] transform = src.window_transform(core_win) 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 # Physical colorbar bounds (real units, not 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 output 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('Easting (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('Northing (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, 'W'), (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"Extent: {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, **({'speed': AVIF_SPEED} if pil_format == 'AVIF' else {})) # 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" {ext.upper()} conversion error: {e}", exc_info=True) return None def _write_subtiles_from(img, tif_file, vis_dir, resolution, subtiles_dir): """Sub-tiles of a tile from its original image (best-effort: on failure, build_index will cut them from the tile image).""" from .index import VIZ_LABELS, _subdivision_k, write_subtiles try: k = _subdivision_k(resolution) stem = Path(tif_file).stem viz_key = next((v for v in sorted(VIZ_LABELS, key=len, reverse=True) if stem.endswith(f"_{v}")), None) if k <= 1 or viz_key is None: return write_subtiles(subtiles_dir, Path(vis_dir).name, viz_key, img, k) except Exception as e: logger.warning(f" Sub-tiles not written ({Path(tif_file).name}): {e}") def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, output_format='avif', subtiles_dir=None): """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). 100 requests lossless, which only WebP honors (Pillow ignores lossless=True in AVIF). output_format: Output format ('webp' or 'avif'). subtiles_dir: pipeline output directory: the map sub-tiles (index_subtiles) are written there from the original image, a single lossy encoding (build_index finds them up to date). Returns: Path to output image file, or None on failure. """ if not tif_file or not tif_file.exists(): return None lossless_gray = any(k in tif_file.stem for k in LOSSLESS_GRAY_KEYWORDS) ext = 'webp' if lossless_gray else ('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 RGB_KEYWORDS) if is_rgb: data = src.read([1, 2, 3]) data = np.moveaxis(data, 0, -1) else: data = src.read(1) # Nodata → NaN: otherwise edge pixels (e.g. -9999, 3.4e38) # pollute the colormap calibration percentiles if src.nodata is not None: data = data.astype(np.float32, copy=False) data[data == src.nodata] = np.nan # Edge stitching: crop to the nominal 1 km tile core_win = _core_tile_window(tif_file, src) if core_win is not None: rows = slice(core_win.row_off, core_win.row_off + core_win.height) cols = slice(core_win.col_off, core_win.col_off + core_win.width) data = data[rows, cols, :] if data.ndim == 3 else data[rows, cols] # 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 from the TIF, or 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 lossless_gray: # Lossless grayscale WebP: exact and 3× lighter than AVIF q100 # (Pillow ignores `lossless=True` in AVIF: lossy q75) img = PILImage.fromarray(rgb_data[:, :, 0]) img.save(str(output_file), format='WEBP', lossless=True) elif quality >= 100: img.save(str(output_file), format=pil_format, lossless=True) else: img.save(str(output_file), format=pil_format, quality=quality, **({'speed': AVIF_SPEED} if pil_format == 'AVIF' else {})) # Map sub-tiles from the original image (after the tile: newer than # it, so build_index does not re-encode them) if subtiles_dir is not None: _write_subtiles_from(img, tif_file, vis_dir, resolution, subtiles_dir) # 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" {ext.upper()} crop conversion error: {e}", exc_info=True) return None def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): """Generate an A3 PDF report for a LiDAR file with all visualizations. Legacy: no longer called by the pipeline (the map's PDF export lives in export_pdf.py). Page 1: Location context (IGN ortho + topo side by side) Pages 2+: Other visualizations (2 per page) Args: basename: Base name for the report file. vis_dir: Directory containing AVIF/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" → Generating A3 PDF report: {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" ✗ No image found for {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: location context 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"Location context - {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" ✓ PDF report done ({time.time()-t0:.1f}s)") return pdf_file except Exception as e: logger.error(f" ✗ PDF error: {e}", exc_info=True) return None